{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cbc3a074-25d9-4a31-be91-d9975aefd062",
   "metadata": {},
   "source": [
    "### 向量自回归模型（VAR）+ AutoEncoder 部分特征"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0844df40-6e6c-4def-82f7-29960a630ec1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from statsmodels.tsa.stattools import adfuller\n",
    "from statsmodels.tsa.api import VAR\n",
    "from statsmodels.tsa.stattools import grangercausalitytests\n",
    "from statsmodels.tsa.vector_ar.vecm import coint_johansen\n",
    "from statsmodels.stats.stattools import durbin_watson\n",
    "\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "import torch\n",
    "\n",
    "from AuEnClass import AutoEncoder\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "pd.set_option('display.float_format', lambda x: '%.4f' % x) \n",
    "np.set_printoptions(precision=5, suppress=True) \n",
    "\"\"\"中文显示问题\"\"\"\n",
    "plt.rcParams['font.family'] = ['sans-serif']\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']\n",
    "plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8bf9f11c-c16c-4139-aae3-7a2025d4e2d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 数据处理\n",
    "data = pd.read_csv('data/train.csv', header=0, encoding='utf-8')\n",
    "df_time_col = data[['时间']]\n",
    "cols = list(data.columns)\n",
    "cols.remove('时间')\n",
    "df = data[cols]\n",
    "# scaler = StandardScaler()\n",
    "# scaler = MinMaxScaler()\n",
    "# df_s = pd.DataFrame(data=scaler.fit_transform(df.values), columns=df.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "527b8df5-a2c5-4c17-8283-54478ce75e47",
   "metadata": {},
   "outputs": [],
   "source": [
    "cols_auto = list(df.columns)\n",
    "cols_auto.remove('总理赔金额')\n",
    "cols_auto.remove('总出险数量')\n",
    "df_auto = df[cols_auto]\n",
    "# scaler_auto = MinMaxScaler()\n",
    "# df_auto = pd.DataFrame(data=scaler_auto.fit_transform(df_auto.values), columns=df_auto.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3f09002c-301e-4471-9bbd-81d80d482178",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch: 001, Training loss= 34858.5977\n",
      "Epoch: 051, Training loss= 34530.2461\n",
      "Epoch: 101, Training loss= 32421.2363\n",
      "Epoch: 151, Training loss= 26793.0156\n",
      "Epoch: 201, Training loss= 19860.1738\n",
      "Epoch: 251, Training loss= 13961.8682\n",
      "Epoch: 301, Training loss= 9584.2646\n",
      "Epoch: 351, Training loss= 6136.9414\n",
      "Epoch: 401, Training loss= 3846.9580\n",
      "Epoch: 451, Training loss= 2834.7830\n",
      "Epoch: 501, Training loss= 2492.7935\n",
      "Epoch: 551, Training loss= 2364.2065\n",
      "Epoch: 601, Training loss= 2312.0051\n",
      "Epoch: 651, Training loss= 2291.1589\n",
      "Epoch: 701, Training loss= 2282.2537\n",
      "Epoch: 751, Training loss= 2277.4680\n",
      "Epoch: 801, Training loss= 2274.0361\n",
      "Epoch: 851, Training loss= 2271.0879\n",
      "Epoch: 901, Training loss= 2268.3792\n",
      "Epoch: 951, Training loss= 2265.8479\n",
      "Epoch: 1001, Training loss= 2263.4783\n",
      "Epoch: 1051, Training loss= 2261.2644\n",
      "Epoch: 1101, Training loss= 2259.2026\n",
      "Epoch: 1151, Training loss= 2257.2871\n",
      "Epoch: 1201, Training loss= 2255.5110\n",
      "Epoch: 1251, Training loss= 2253.8674\n",
      "Epoch: 1301, Training loss= 2252.3477\n",
      "Epoch: 1351, Training loss= 2250.9434\n",
      "Epoch: 1401, Training loss= 2249.6453\n",
      "Epoch: 1451, Training loss= 2248.4453\n",
      "Epoch: 1501, Training loss= 2247.3342\n",
      "Epoch: 1551, Training loss= 2246.3030\n",
      "Epoch: 1601, Training loss= 2245.3430\n",
      "Epoch: 1651, Training loss= 2244.4470\n",
      "Epoch: 1701, Training loss= 2243.6064\n",
      "Epoch: 1751, Training loss= 2242.8145\n",
      "Epoch: 1801, Training loss= 2242.0632\n",
      "Epoch: 1851, Training loss= 2241.3479\n",
      "Epoch: 1901, Training loss= 2240.6628\n",
      "Epoch: 1951, Training loss= 2240.0022\n",
      "Epoch: 2001, Training loss= 2239.3621\n",
      "Epoch: 2051, Training loss= 2238.7380\n",
      "Epoch: 2101, Training loss= 2238.1272\n",
      "Epoch: 2151, Training loss= 2237.5259\n",
      "Epoch: 2201, Training loss= 2236.9319\n",
      "Epoch: 2251, Training loss= 2236.3428\n",
      "Epoch: 2301, Training loss= 2235.7566\n",
      "Epoch: 2351, Training loss= 2235.1719\n",
      "Epoch: 2401, Training loss= 2234.5876\n",
      "Epoch: 2451, Training loss= 2234.0020\n",
      "Epoch: 2501, Training loss= 2233.4150\n",
      "Epoch: 2551, Training loss= 2232.8259\n",
      "Epoch: 2601, Training loss= 2232.2339\n",
      "Epoch: 2651, Training loss= 2231.6384\n",
      "Epoch: 2701, Training loss= 2231.0403\n",
      "Epoch: 2751, Training loss= 2230.4382\n",
      "Epoch: 2801, Training loss= 2229.8323\n",
      "Epoch: 2851, Training loss= 2229.2227\n",
      "Epoch: 2901, Training loss= 2228.6096\n",
      "Epoch: 2951, Training loss= 2227.9924\n",
      "Epoch: 3001, Training loss= 2227.3713\n",
      "Epoch: 3051, Training loss= 2226.7471\n",
      "Epoch: 3101, Training loss= 2226.1189\n",
      "Epoch: 3151, Training loss= 2225.4878\n",
      "Epoch: 3201, Training loss= 2224.8533\n",
      "Epoch: 3251, Training loss= 2224.2151\n",
      "Epoch: 3301, Training loss= 2223.5752\n",
      "Epoch: 3351, Training loss= 2222.9316\n",
      "Epoch: 3401, Training loss= 2222.2859\n",
      "Epoch: 3451, Training loss= 2221.6379\n",
      "Epoch: 3501, Training loss= 2220.9873\n",
      "Epoch: 3551, Training loss= 2220.3350\n",
      "Epoch: 3601, Training loss= 2219.6799\n",
      "Epoch: 3651, Training loss= 2219.0215\n",
      "Epoch: 3701, Training loss= 2218.3618\n",
      "Epoch: 3751, Training loss= 2217.7007\n",
      "Epoch: 3801, Training loss= 2217.0378\n",
      "Epoch: 3851, Training loss= 2216.3743\n",
      "Epoch: 3901, Training loss= 2215.7095\n",
      "Epoch: 3951, Training loss= 2215.0437\n"
     ]
    }
   ],
   "source": [
    "# AucoEncoder\n",
    "# 更改网络参数\n",
    "# input_dim = 13\n",
    "# hide_dim1 = 6\n",
    "# code_dim = 2\n",
    "\n",
    "n = df_auto.values.shape[0]\n",
    "autoencoder = AutoEncoder()\n",
    "autoencoder.train(df_auto.values, epochs=4000, batch_size=n)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f55ded75-6f6f-4ea6-b11b-34e79fe0f967",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 利用训练好的自编码器重构数据\n",
    "# 数据类型转换\n",
    "temp_data = torch.FloatTensor(np.array(df_auto.values))\n",
    "encodedData, _ = autoencoder(temp_data)\n",
    "encodedData = encodedData.double()\n",
    "encodedData = encodedData.detach().numpy()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e5a46bdf-cef7-4d0e-b3c6-6b6373c31996",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 总理赔金额，总出险数量 进行归一化\n",
    "df_res = df[['总理赔金额', '总出险数量']]\n",
    "# scaler_res = MinMaxScaler()\n",
    "# df_res = pd.DataFrame(data=scaler_res.fit_transform(df_res.values), columns=df_res.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "79cc7fdb-ad05-4594-af4a-3d18155f922b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(data=np.concatenate((encodedData, df_res.values), axis=1), columns=['V1', 'V2', '总理赔金额', '总出险数量'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "edd88918-15e7-4cfc-9428-84832791ee22",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAAGoCAYAAACg4OZuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd5xU9dXH8c+hKSUiyIqxIKLYuwRFQbE3TCx51MQaY/AxJs9jNEWj6bZoYrpGLLFETbDG3kUQRYXHAtijYEVpUsRCOc8fZyYMy+zuzO69M3d2vu/Xa17M3pm5c/bu7uXMuef3+5m7IyIiIiLS3nWodgAiIiIiIpWgxFdERERE6oISXxERERGpC0p8RURERKQuKPEVERERkbqgxFdERERE6oISX2lXzOxJM/uvgq/PN7NLcvePMrO/VS86EZH2ralzsJmdYWbzzOwjMzu9mjFKfVPiK+3Nw8Cwgq+HAg+b2QjgL4BVJSoRkfpQ7Bz8LvA1YDtgZ+BnZrZeFWITUeIr7c4j5E66ZrYKsAMwBvgG8KvqhSUiUheKnYMfB77u7m+4+4vANGDtqkUodU2Jr7Q344GBZtYTGAy87O6zga8Cs6samYhI+1fsHPyYu08FMLO1gfWAKVWMUepYp2oHIJIkd//MzCYQl9O2Iy674e5upi4HEZE0NXUOLnAucJm7f1zx4ERQ4ivtU/5S23bAH6oci4hIvSl6Djaz/Yme322rFJeIWh2kXXoY2AP4EjCuyrGIiNSblc7BucFsVwBHqdor1aTEV9qjicCmwEs6wYqIVNwK5+DcILc7gIvc/enqhib1TomvtDvuvhR4jLjcJiIiFVTkHLwf0d5whpnNyN0Oq1qAUtfM3asdg4iIiIhI6lTxFREREZG6oMRXREREROqCEl8RERERqQtKfEVERESkLlRlAYs+ffp4//79q/HWIiIVM2nSpFnu3lDtOArp/Csi9aCp829VEt/+/fszceLEary1iEjFmNn0asfQmM6/IlIPmjr/qtVBREREROqCEl8RERERqQtKfEVERESkLiSa+JrZlWb2pJmdneR+RURERETaKrHE18wOBTq6+xBggJkNTGrfIiIiIiJtlWTFdzgwOnf/AWBo4YNmNtLMJprZxJkzZyb4ttIuzZgB06ZVOwoRqYRJk+CTT6odhYjUgSQT3+7Au7n7c4C+hQ+6+yh3H+TugxoaMjWtpWTR978PJ51U7Sjqy+jRMHdutaOQenTKKTBuXLWjEJE6kGTiuxDomrvfI+F9Sz1ZsADuuguefhrefbfl50vbffopnHBCJL8ilfbhhzB/frWjEJE6kGRyOonl7Q3bANMS3LfUk5tvht12g8MPh+uuq3Y09eGxx2DxYrjzzmpHUtfMrLeZ7W1mfaodS0XNmhUfeEVEUpZk4ns7cIyZXQwcDtyd4L6lnlxzDRx3XNyuvhrcqx1R+3fnnfC978HYsbBoUbWjqUtm1gu4CxgMPGpmDWb2lpmNyd22yj1vpdlzanpGnc8+i6RXFV8RqYDEEl93n08McJsA7O7u85Lat9SRadNgyhQ48EAYMgSWLYuWB0mPe7SWHHMM7LADPPJItSOqV1sDp7n7ucD9wAnAje4+PHebXGz2nJqfUSc/2FkVXxGpgE5J7szd57J8Zof6MW8eTJ8Os2fHJbvC2+zZ8NFH8LOfwZe+VO1Is++66+CII2CVVeLrfNV3xx2rGla7NnUqmMHmm8OIEZEEjxhR7ajqjrs/BmBmuxJV35uBEWa2OzAZOInis+dsV2TbaxULvK3yia8qviJSAYkmvu3e3Lnw4ovLb1Onxr8ffQT9+0NDA/TpE7c11oANNohk98034TvfgQkTIsGQ4tzh2mvh+uuXbzvmGNhuO/jd72DVVasXW3t2111w0EHxuzliBOy5Z/ws9LtacWZmwBHAXOBZYC93f9/MrgUOYOXZc7ZvYlvj/Y4ERgL069cvzW+hfKr4ikgFKfEtxezZkXzNnRtVsc03hy22gL33jvv9+kGHZrpGli2D226DW2+Fww6rXNy15oknoFOnFSvj/frB9tvDHXfEYDdJ3l13wU9/Gvc33hi6doXnnovfeakod3fgFDP7FbC2u+fn+JoIDKT47Dktzqjj7qOAUQCDBg3KVtP8rFnxIUsVXxGpACW+pXjvPejePdoZWlMF69ABLrgAvvtd+PKXoXPn5GNsD669Fo49duVjnG93UOKbvFmzYPLkmEUD4tgfdFAkw0p8K8rMfgS87+7XAqsDfzWzl4EpwMHAecBMopVhAjF7zivAO0W21Y6ZM2G99VTxFZGK0Fy7pViwAFZfvW2XfvfeO07uV12VXFztySefwE03RWtDY4ccAk8+Ce+/X/m42rv77ovWhnxPNSzv85VKG0XMjDMW6AjsClwHPAc86e4PUXz2nNqeUWfmTBgwQImviFSEEt9SzJ8PX/hC2/ZhFlXfX/wCPv44mbjakzvuiBkF1l135ce6d4dDD4W//73ycbV3d9658kC2oUPh1Vdj2WipGHef6+57u/uu7v5td5/s7lu7+1buflbuOSvNnlPzM+rkE1+1OohIBSjxLcWCBbDaam3fz6BBsOuu8Ic/tH1f7U1+7t6mHH98PEdz+iZn8WJ44AE44IAVt3fpAvvsA/fcU524pFm5BHm0u89oblvNUMVXRCpIiW8pkqj45p1zDlx8cfRWSpgxI1oZDjmk6ecMHRrtEJMmVS6u9u7xx2HgQFhrrZUfU7uDVMqsWbDhhqr4ikhFKPEtRVIVX4CNNoIjj4Tzzktmf+3B9dfDwQdHS0NTzJYPcpNkNDdf7/77w8MPx6paImlSxVdEKkiJbymSrPgC/OQncdl+2rTk9lmr3Ftuc8g79lj4xz+UjCWlucS3Tx/YcksYM6aiIUkdmjkz5kH/9FNYsqTa0YhIO1c7ie+SJfDoo9V57yQrvgB9+8aCFvm5U+vZ88/HB4tdd235uf37w1Zb6RJ8El59FRYubH7Ksvy0ZiJpWbo0FgBaYw3o0SN+J0VEUlQ7ia8ZHHhgdU6MSVd8AU4/PQYWPf98svutNddcE1OYNbcASKH8IDdpm3y1t7kp+vJ9vhpQKGmZPTumiuzYMc6x6vMVkZTVTuLbsWOsKvXyy5V/76QrvhD7O+ssOPPMZPdbSxYvhhtuiBaGUh12GIwdCx98kF5c9aC5Noe8LbaIpPfFFysTk9SfWbNiqXeIc6L6fEUkZbWT+EIsDzx1auXfN42KL8BJJ0UiX699lPfdF4P9Bg4s/TU9esRAuOuvTy+u9u6jj2DixFi4ojn5VdzuvLMycUn9mTlzeeKriq+IVECbE18z62lm95rZA2Z2m5l1SSKworbYojqJbxoVX4j5Us89F37wg/qsdFx7bWmD2hrLz+6gS/Ct88ADMGwYdOvW8nM1rZmkqXHiW4/nQRGpqCQqvkcBF7v7PsAMYL8E9lncFltU57JrWhVfgCOOgK23jjaOP/2pfmYsmDMHHnwQDj+8/Nfutlv8TJ57Lvm46kEpbQ55u+0Gkydr3mlJx8yZMYMIqNVBRCqizYmvu1/i7g/mvmwAPmzrPptUrVaHtCq+EIO6rrwyLvvfdx9stlkszbtsWTrvlxXXXRcrhq2+evmv7dAh5kK+6abk42rvli6Fe++NgaKlWHVV2GOP+N0USZpaHUSkwspOfM3sMjMbU3D7aW77EKCXu09o4nUjzWyimU2cOXNm66LdcMMY1PTxx617fWulWfHN22YbuPvuuIR/ySUxzdTdd7fPy/nucNll0ePcWl/9Ktx8czrHp1rziX7+OVx4Ibz2Wnrv8dRTsPba0K9f6a9Rn6+kRYPbRKTCyk583f0kdx9ecPulmfUG/gSc0MzrRrn7IHcf1JA/0ZUrP7PDSy+17vWttWBB+olv3q67wvjx8Mtfwg9/GJeax4xJLwFevDid/Tbn8cfj+yll7t6m7LBDJIpTpiQXlzvceGN8wNp778r+J/zKK7DzzvDrX8ciHa3x6qtRyf3nP+PYFHPnnaW3OeQdcED0BVfjd0XaN1V8RaTCkhjc1gW4CTjT3ae3PaQWVLrdYfHiSCJKGQiUFDP4ylfghRfghBPg29+OJT3POiuZHudPP41ZEYYOjVkSLr+87fssx2WXwciRzc8h2xKzmNrs5puTiWnq1Likf+GFkXhuskl8nXZvq3sc/6FD4ZvfhKuuig8GrXHXXbBoURzf9dePFQLffnvl55Sb+K61Vsy80dq4pGRm1tvM9jazPtWOpSLU4ysiFZbE4LZvAtsDZ+VaH45IYJ9Nq/QAt3y1ty1JWmt17BgLNkydCrfeGgn4PvtEG8RvfgPvvlve/l5/PWaQ6NcvFoE4/XT4v/+D3/8eTjwxEuK0zZ4dyVdrZnNoLN/u0BYLFsD3vw/Dh0ci/cwzMePBpZfCXntFxb3c41yqWbPg0EOjtWXsWDj55EiAJ0xoXavF2LHRPvLII3GbNw+23Tbe46GH4M03o1Vo8ODy961V3FJnZr2Au4DBwKNm1mBmV5rZk2Z2dsHzStpWE1TxFZEKS2Jw26Xu3qug9eGfSQTWpEpXfNMc2FYqs0h2L7oIpk+Hiy+O+X+32gp23x1+9jP44x9jMYj774dJk2DatIh98eJImvfZB4YMif098URcuj7kkPgg8dRT8R/OsGHw1lvpfi/XXBNJVO/ebd/XjjvGnLStaX3JtzVsumnMMDF1aiwj3alTPG4G558fi2sMGxYfGpL04IORlG64YSS6m20W29dYA9ZdN6r95Vi2DMaNW94+stlm8TsxfTrsuy+cdlr8Dh1wQHygKteIEXDbbVFRrrS3345EfunSyr93ZW0NnObu5wL3A3sAHd19CDDAzAaa2aGlbKvad1CuwsRXFV8RqYBO1Q6gbJWey7cSA9vK0bFjJLu77w5//jPcc08kSa++Ck8+GVXE2bOX//vpp5Hwnnwy3HFHjNJvrEeP6Au9+OKoBl5/fcuLG7RGflDbVVcls78OHaJKe8stcHYZha4ZM2JWiHnzomKc/0BQzI9+BL16ReX33ntj6rm2+Owz+PGP43hffXVUlRsbNiyS2O23L32/U6fGh4m1115xe48eUQUeOTJ+Pxo/Xqptt4Vddokq++23xxzUaXrrrfjZ3HRT/G6vuy588klcpTj2WOjaNZn3eewx+MtfIulad93lt3XWiX9XX71iV3vc/TEAM9uVqPr2BkbnHn4AGApsV+K2FUZImtlIYCRAv3IGNqbJfcXBbar4ikgF1NbKbRC9rjNmVG5mhyxUfJuy6qpxGfvnP48k+MYbo5L4f/8XicPHH0ey8PjjcNRRxZPePLNIKm64AY4+Onpdkx5QN2YMdO4cg7iS0pp2h3PPjR7eiRObT3rzRo6MDwV77x3JY2stWhSV09deg+efL570QrQ7jBtX3r7Hjm1+sKBZHPf+/cvbb+Hr//a3SHiPPjqd6uv06fDb38JOO0XSP3VqXM14//2Ys/mKK6LdYoMN4JxzolLfWlOnxpWH44+PXu7Bg6Nq/uST8LvfxQej/v3jg8M3vpHUd9giMzPgCGAu4EC+z2YO0BfoXuK2FSQyuDhp8+fHOWmVVeJrVXxFpAJqr+LbqdPymR0GDUr//bJW8S1XuZW5PfaAp5+OSurTT0eyk9T3n5/CLMkK2i67xAeh11+P5Y9b8uGHUdF+8cXyLvkfcQT07BmDDq+/PpLgcixaBF/+MvTtG+0enZr50xs2LPqO3Us/VmPHRhtDmjp1ioF/I0bEh4HLL4+qe1vMmRP7vO66+EBw8MHwi1/E72Hnzis+d9dd4zZ1avS4b7QRHHMMfO97pSf0770XyfS//gVnnhkfmvKJVzELFlR0+kR3d+AUM/sV8FUgP/K0B1GoWAh0LWFb9hUObANVfEWkImrjBNlYJQe4Zbnim5b11otEqnfv6EH99rfj67YsqvHhh7EIwjHHJBcnRPJ6yCHR7lCKP/4xkti11ir/vfbbL/qljz4a/vCH0o/HokVRXfziF2OZ5uaSXojBh507l95X7N5yxTcpq64arQ4vvRRXCFpzVWDx4mi7OeywqN6OHRszULz/flR199135aS30BZbxAeyyZMjnh12iET5f/8XRo2KHvZ581Z8zfz50Q6z1Vbxe/3KK5EwN5f0QiRjrfldaQUz+5GZHZv7cnXgAqJtAWAbYBowqcRt2VfY3wtaslhEKqJ2E99K9fnWesW3tVZdNZKICRMiEf7udyMhO/30mPmg3ITn6qsjQW3NSm0tKbXdYf58+OtfY2aL1ho6NOZZvvHGqLC+917zz//446iQrrNOHINSqsxmUfUtdfqw11+PRLG1bQzl6tEjFld59NGozpbCHZ59Fk49NXpnL7wwPkhMnx4V3wMOaD7ZLWaddWLe4zffjDmv+/WLVoVTT43H1lsv3uPb346rRG+/HTH8+tfRt509o4BjzGws0BG4Pff1xcDhwN1lbMu+wv5eUKuDiFRE7bU6QMzscOWVlXmveqz4FhowIC4Jn3lmVNn/8Q/4+tcjkTnyyEiI+67UUriiZcsiif7739OJcbfdIvmZNq355O+yy2J2iwED2vZ+G20USek550Qv6qWXRlLfWD7pXX/9+H0tp7UiP8CtlP7Sxx6Lam8lp9zr1StmBhk2LFpAvve9lZ/jHr25t94aFflFi2Jg2vjxpbWllGq11SLB3W+/5duWLYukesqUqO7ed18M0Mswd58LrNBDY2bDc9sudPd55WzLvGIVX7U6iEjKarfiW6lWh3qt+Baz+eaxotyrr8asBLNnx4CplpbYfeQR6N49ph9LQ6dO0Rt6661NP+fTT2PQ0hlnJPeeP/95TPH1gx/E4hMLFy5//OOPYxW1/v3LT3qhvAFulWpzaGzNNWN+4D/8IVoUIBLO8ePjysCAAfBf/xUzWVx1FbzxRvz+JJn0NqVDh2ijOOig6JfOeNLbFHef6+6j3X1Gudsyr3Hi2717/J22/2nrRKSKajPxHTAg+gErMeik3iu+xZhFX+Wll0YleLfdYu7gpqQxqK2xltodrr025rFt63RkjQ0ZEpfPIZKrCRMiAT7ggPg9bU3SC/HhbvbsGLjXkmolvhDtBA8+GAPGvva1aDE4+eT4sPivf8WHogsvjJka2joQTtqXxoPbOnSI5FftDiKSotr8nyg/s8PLL6f/Xqr4Nu/EEyMB3n//SIAamzEjqoJHHZVuHHvsEb8P77yz8mNLlkTydeaZ6bz3F74QCe6FF0bledCgWOL3iitan+x16BDV9Jb6fKdPjxaCTTZp3fskYeDA+BnvuGMk4S+8ENXwrbeuzoqHUhsaV3xBfb4ikrraTHyhcgPcVPFt2Ve+Ej2cRx8dPcCFrroqqrE9e6YbQ5cucVn7tttWfuyWW2Jk/tChKz+WpEMPjTmUTz01eprbWuEsZYBbfrW2aieYm20W3/fA2lk0TKqs8eA2UJ+viKSudhPfSi1drIpvaYYNi6rf978fU4ZB9Htefnm0OVRCsXYH91h6OK1qb2Nrrw3//d/JXNbPD3BrzmOPRauJSK1RxVdEqqB2E99KDXBTxbd0W20VFcq//CWW5b3//pgztRILjUAsKvH88yv2xd5/fyTgaS/ukIYddogZCZqrgFWzv1ekLYolvqr4ikjKajvxVcU3e/r3j1H9Dz8cg50qVe2FmHv4gANigYW888+PmRyq3QrQGqusEsnvhAnFH58xIxYG2XLLysYlkoTGg9tAFV8RSV3tJr6VmtlBFd/y9ekTie93vhNz/lZSYbvDE0/EogWHH17ZGJLU3LRm48bF462ZNUKkmj75JFbwa1xUUMVXRFKWWOJrZn3N7Nmk9teiSs3soIpv6/ToEQs89OhR2ffdb79YWW7mTLjggphjt6UlgrOsuQFuanOQWpUf2Nb4SowqviKSsiQrvr8Buia4v5ZVYoCbKr61pVu3WJ3tvPMiAT7++GpH1DZDhsT38fnnKz+mgW1Sq4r190IUGZT4ikiKEkl8zWwP4GOgsqsGpT3AbfHiuHWtbD4vbfTVr8Lvfw//+7+1/7Pr2TOmCGu8QMicObFE83bbVSUskTYp1t8LanUQkdSVnfia2WVmNqbg9lPgJ0Cza8Ga2Ugzm2hmE2fOnNnaeFeU9gC3BQviRFyLA6Pq2YEHwvDhsYJYe1Cs3eHxx2M1tM6dqxOTSFs0VfFVq4OIpKzsxNfdT3L34flbbvMl7v5RC68b5e6D3H1QQ7ETXmuk3eqg/t7a1KMHPPpo+otmVEqxAW7q75Va1lyrgyq+IpKiJFod9gJOMbMxwLZmdkUC+yzNhhumO7OD+nslC4YOjSnili1bvk2Jb7tjZj3N7F4ze8DMbjOzLmb2VsHVta1yz7vSzJ40s7MLXrvStkwrtmobqOIrIqlrc+Lr7rsWVH+fc/cT2x5WidKe2UEVX8mCtdeG1VeHl16KrxcsiN72wYOrG5ck7SjgYnffhxgvcQZwY8EVtslmdijQ0d2HAAPMbGCxbdX7Fkqkiq+IVEmi8/gWtD5UzuabpzfATRVfyYrC5YufeCIWtlh11erGJIly90vc/cHclw3AEmCEmT2dq+h2AoYDo3PPeQAY2sS2FaQyxqItmhrcpoqviKSsdhewyEtzgJsqvpIVhQPc1ObQrpnZEKAX8CCwl7sPBjoDBwDdgXdzT50D9G1i2wpSGWPRFqr4ikiV1H7im+YAN1V8JSsKB7gp8W23zKw38CfgBOAFd38/99BEYCCwkOXzpfcgzuHFtmWbZnUQkSrJ/gmyJWnO5Zufzkyk2jbeOJZ5ffVVePbZWNhC2hUz6wLcBJzp7tOB68xsGzPrCBwMPA9MYnkrwzbAtCa2ZVtTg9u0gIWIpKyG13LN2XBDeO89WLQoVu1K0vz5qvhKNphF1fe3v4Utt6z8UtBSCd8EtgfOMrOzgEeB6wAD7nD3h8xsNWCcma0N7A/sBHiRbdm1ZAnMmwe9eq38WPfu8QFv6VLo2LHysYlIu1f7Fd9OnWJlqzRmdlDFV7Jk2DC4+mq1ObRT7n6pu/cqmMXhF+6+tbtv5e5n5Z4znxjMNgHY3d3nFdtWre+hJLNnR9JbLLHt0CGS34ULKx+XiNSF2k98Ib0Bbqr4SpYMHQqff67Et865+1x3H+3uM5rblllN9ffmaYCbiKSofSS+aQ1wU8VXsmS77eJ3fZddqh2JSOu1lPhqgJuIpKh9JL5pDXDTdGaSJZ06xQe8Yr2RIrWiqYFtear4ikiK2k/im1bFV60OIiLJUcVXRKqofSS+hTM7JEkVXxGRZDW1alueKr4ikqL2kfimNbODKr4iIslSxVdEqqh9JL6QTruDKr4iIslSj6+IVFHtL2CRt/nmKw5wW7YMPvssJkP/5BP49FPo1w86dy59n6r4iogkSxVfEami9lPx3XZb+N3vYsT7qqvG5Oi9ekX/75e+BNtvD3/9a+n7W7w4bl27pheziEi9KaXHV4mviKQksYqvmV0C3Ovudya1z7KMGAFvvgmrrBLJ6iqrxCpAeeedFwPgSpWfw9cs+VhFROpVKRXft9+uXDwiUlcSSXzNbBiwVtWS3ggCvvjFph9vaIB//7v0/am/V0QkWe7R46uKr4hUSZtbHcysM3A5MM3MvtL2kFKy5ppRaSiV+ntFRJI1bx506xZX5Jqy2moa3CYiqSm74mtmlwGbFGx6FHgRuBD4rpn1c/c/FXndSGAkQL9+/VoXbVs0NMCHH5b+fFV8RUSS1VKbA6jiKyKpKrvi6+4nufvw/A1oAEa5+wzg78DuTbxulLsPcvdBDS2d+NKgiq+IZJiZ9TSze83sATO7zcy6mNmVZvakmZ1d8LyStmVSSwPbQNOZiUiqkpjV4XVgQO7+IGB6AvtMniq+IpJtRwEXu/s+wAzgSKCjuw8BBpjZQDM7tJRtVfsOWlJKxVfTmYlIipIY3HYlcJWZHQl0Br6awD6Tt9pq8PnnMadvKVOUqeIrIhXk7pcUfNkAHA38Pvf1A8BQYDtgdAnbXks73lYptdVBFV8RSUmbK77uvsDd/8vdd3X3Ie7+bhKBJc4sTriltjuo4isiVWBmQ4BewNtA/nw6B+gLdC9xW+N9jjSziWY2cWY5LV9Ja2nVNlDFV0RS1X4WsChFOX2+qviKSIWZWW/gT8AJwEIgf3mqB3G+LnXbCqo+xiKvlIpv9+6waBEsXVqZmESkrtRX4ltOn68qviJSQWbWBbgJONPdpwOTiLYFgG2AaWVsy6ZSBrd16AA9esDChZWJSUTqSmIrt9WEciu+G22UbjwiIst9E9geOMvMzgL+BhxjZmsD+wM7AQ6MK2FbNpVS8YXlU5r17Jl+TCJSV1TxbUp+yWIRkQpw90vdvVfBdJHXAMOBCcDu7j7P3eeXsq0630EJSk18tYiFiKREFd+mzJ+vHl8RqSp3n8vyGRvK2pZJpQxuAy1iISKpqa+K75prquIrIlItqviKSJXVV+Jb7nRmqviKiCQjP1ND9+4tP1cVXxFJSX0lvqr4iohUR77aa9byc7WIhYikpL4SX1V8RUSqo9Q2B9AiFiKSmvpKfFXxFRGpjlIHtoEqviKSmvpKfLt3h2XL4OOPm3/e4sVx69q1+eeJiEhpSlm8Ik8VXxFJSX0lvmalTWmWr/aW0osmIiItK6fVQRVfEUlJfSW+UFqfr/p7RUSSpR5fEcmA+kt8S+nzVX+viEiyyq34KvEVkRS0OfE1s15mdo+ZTTSzy5IIKlWq+IqIVF45g9u0gIWIpCSJiu8xwPXuPgj4gpkNSmCf6VHFV0Sk8soZ3KaKr4ikJInEdzawpZmtDqwHvJ3APtOjiq+ISOWV2+Oriq+IpKBTuS/ItTNsUrDpUWB94H+Al4A5TbxuJDASoF+/fmUHmpg114SpU5t/jiq+IiLJUo+viGRA2RVfdz/J3Yfnb0TS+9/u/kvgZeAbTbxulLsPcvdBDaWe/NLQ0NByq8P8+Up8RaTizKyvmY3L3e9kZm+Z2Zjcbavc9ivN7EkzO7vgdStty5TFi2HhQujVq7Tnq+IrIilJotWhF7CVmXUEdgQ8gX2mp9R5fNXqICIVZGa9gGuA7rlNWwM3FhQaJpvZoUBHdx8CDDCzgcW2Vec7aMbs2dC7N3Qo8V4gt4MAACAASURBVL+c7t3hk09g6dJ04xKRupNE4ns+MAqYB/QGbkxgn+lRxVdEsmkpcASQL3XuBIwws6dzFd1OwHBgdO7xB4ChTWxbgZmNzM28M3FmSx/801DOwDaIBLlbt6gSi4gkqM2Jr7s/7e5buHsPd9/b3bN9psoPbvNmCtOq+IpIhbn7fHefV7DpGWAvdx8MdAYOIKrB7+YenwP0bWJb431Xt9WsnP7ePC1iISIpKHtwW83r3j2qCQsXNl3VVcVXRKrvBXf/LHd/IjAQWAh0zW3rQRQvim2rrI8+ggsugA02gM02g003jUQ3v+x7axJfLVssIimov8QXlvf5NpXcquIrItV3nZmdC0wBDgbOA2YSrQwTgG2AV4B3imyrrHHj4F//gp13hmuugZdeigLDZpvFbeZMWHvt8vapiq+IpKA+E998n++AAcUfV8VXRKrvl8ANgAF3uPtDZrYaMM7M1gb2J/qAvci2ypoyBUaMgIsuiq/d4YMPIgHO3w46qLx9akozEUlBfSa+Lc3soIqviFRJbppI3H0KMbND4WPzzWw4sDdwYb4nuNi2ipo8Gfbdd/nXZrDWWnHbfffW7VNTmknWTJkCW2yxvIVHalLle8GyoKWZHbSAhYhklLvPdffR7j6juW0VNWUKbLVVsvtUxVeyZNkyGDoUxo+vdiTSRvWZ+LZU8dWSxSIipVm8GF57LXp5k6SKr2TJ1Kkwbx48/ni1I5E2qs/EVxVfEZFkvPYarLcedO3a8nPLoYqvZMn48ZE7jBtX7Uikjeoz8W2u4rt4cdySPomLiLRHkyfDllsmv19VfCVLxo+H//1feOIJrShY4+oz8W2u4puv9qp5XUSkZWn094IqvpIt48fDIYdE/jB1arWjkTaoz8S3uYqv+ntFREo3ZUo6FV8tYCFZ8f77sUjLppvCsGFqd6hx9Zn4llLxFRGRlqWV+GoBC8mKJ56AIUNiUZahQzXArcbVb+I7c2ZMst6YKr4iIqX5+GN4913YaKPk961WB8mK8eNhl13ifr7iWyx/kJpQn4lv167QpUvxy2iq+IqIlOall2DjjaFz5+T3rcFtkhWFie+GG8KSJTB9enVjklarz8QXmu7zVcVXRKQ0abU5gCq+kg2LFsXv+Ze+FF+bqc+3xtVv4ttUn68qviIipUkz8VXFV7LgmWfid7xbt+Xbhg1Tn28NKzvxNbO+Zjau0bYrzexJMzs7udBSpoqviEjbqOIrxbz1VvvpgS1sc8gbOlQV3xpWVuJrZr2Aa4DuBdsOBTq6+xBggJkNTDbElKjiKyLSNpMnpzOHL0CPHnGZedmydPYv6Zg3DwYOhIcfrnYkySiW+G69dQzqnDWrOjFJm5Rb8V0KHAEUXn8aDozO3X8AGFrshWY20swmmtnEmU3NoVtJqviKSMY0vqJW7GpaqdtSN2dOFAr69Utn/x06xOXlhQvT2b+k4667YgD5b35T7UjabtkyePJJ2HnnFbd36gQ77RRJsdScZhNfM7vMzMbkb8Cp7j6v0dO6A+/m7s8B+hbbl7uPcvdB7j6ooaGhrXG3nSq+IpIhja+oFbuaVuq2igQ8dSpssUW6q1yqz7f23HorXHghvPBC3GrZyy/D6qvDF7+48mPq861ZnZp70N1PKmEfC4Guufs9qJUBc2uuCZMmrbxdFV8RqY78FbV/5b4ezspX07YrcdtrhTs2s5HASIB+SVVo0+zvzVOfb235+GN46CG4/PK4IvDb38I111Q7qtYr1uaQN3QonHFGZeORRCSRpE5ieXvDNsC0BPaZPlV8RSRD3H1+oytqxa6mlbqt8b6Tv+KWZn9vniq+reMON99c+TaR++6DwYOhd2846SS48054553KxpCk5hLfwYPjb2DRosrGJG2WROJ7O3CMmV0MHA7cncA+06ceXxHJtmJX00rdlj5VfLNp2TI4/XQ4/HC46qrKvvctt8Bhh8X9Xr3guOPgj3+sbAxJai7x7dYtBrk99VRlY5I2a9UJ0t2HF9yfT1ySmwDsXqQHOJtU8RWRbCt2Na3UbelyV+KbRUuWwDe/CRMmRK/tqFGVm1bss8/g3nvh4IOXbzv1VLjyytqs2n/wQRTHttii6eeoz7cmNdvjWyp3n8vyHrPa0NAQU5G4rzg4QxVfEcmG24FxZrY2sD+wE+AlbkvX++/HyPY110z3fdTqULrPPoOvfz0+KDz4YFQkf/Sj4rMSpOGhh6L1Za21lm9bf33Yd9/o+T399PRjWLoUnn46Zlxo66DLJ56AIUNidpGmDB0Kf/5z295HKq42BqKlYZVVYsqVjz5acbsqviJSRfkrasWuppW6LfUgJ09Ov9oLqviWauFCGDEikrQ774Tu3SPxGzkSLrusMjHccgsceujK27//ffj972Hx4nTf/+23Ya+9YI894Je/bPv+mmtzyNtll6iuL1nS9vdrq7SPbztSv4kvFO/zVcVXRDLC3ee6+2h3n1HutlRNmZL+wDZQxbcUc+ZEwrf++vCPf0RRJ++44+Bf/4K5c9ONYfFiuOOO4onv9tvDxhvDP/+Z3vvfdBPssAPssw+88UYch/POa9s+n3ii5cR3jTVgvfXg+efb9l5t9dOfxoed/faDK67QwhotqO/Et3Gf7+efxye3VVetXkwiIllXif5eUMW3Je+/D7vtFpfcL78cOnZc8fE+fWD//eHvf083jrFjYcCAphcz+f734aKLku83XrAAjj8ezjoL7r4bzjwz5tx95BG4+urWL6Lx6aeRzA4e3PJzq93n++c/R6L/yivR3/3gg7DhhvFh6K9/jV5lWUF9J76NK74LFkSFIc0J2UVEal2lEl9VfJv2wguRdH3ta5FUNvX/1siR6Q9yK5zNoZj99ovZJh56KLn3nDABtt0WOneG//s/+NKXlj+WT34vvRT+8Ify9z1xImy2WVRRWzJ0KIwb1/Lz0nDTTXD++XD//bDBBvBf/xWV9fffh1NOibg23RSGD4+Wl3m1MfdA2pT4FlZ81d8rItK8ZcvgpZeaH+2eFFV8V/TSS9G/utVWUcn90Y/gxz9uvlgzfHgMfJswIZ2Yli2D224r3uaQZxaD2y66qO3vt2RJHIOvfCX2d/nl0KPHys9bd91Ifn/3O7jkkvLeo5T+3rx8xbdSs2fkPfJIJLd33x1Jb6Fu3eCQQ+D66yMJPu20+NCx/vpw9NHw8MPxc0vLvHlxTC69FL7znVjEZOnS9N6vTPWd+DY0FK/4iohIcW+8EZfQK3GuVMUXXnwRfvGLqLDvvXf09P71rzGY61vfavn1+UFuo0alE98TT8T/pQNbWCn761+PZa7b0g+7cGH08Y4bF1Xe5pJtiETvkUfgggui97VU5SS+/fpF1fn110vff1s9+ywceSSMHh1V7+asuip8+ctRHX799WjfOP30aE35+c9h2rS2xTJ7NtxwQ7SZjBgRx3yddeI9Jk2KpPyKK2LO49tuq/wHhCLqO/FtXPGdP18VXxGR5lSqzQHqu+J7331RVd9335h9aNQoeOutmCFhl12an2arseOOg9tvX3kWoyTcemvzbQ55XbrA//xP6/tuP/44EqsNNohL++usU9rrBgyICufPf17a8snupQ1syzOrbJ/vG2/Ecbj00qjml6NPn/gZPPdcJKFz5sCgQbDnnpG8fvpp6ft6552oJA8cGEl1167RY/zII5FLPfVUJLynnx494BddFB/gdtopnlNFiczjW7MaGla8/KNWBxGR5lUy8a3Hiu/SpfCrX8Ul/KuuiipvOUluMQ0NkUBff31cHk+KeyS+d95Z2vNPOikS0bffjtkQSrVo0fKk9/LLyz8eAwfGoK8994xE9dhjm37uK69E60SpiTUs7/P9xjfKi6tcH3wQP8ezzy7tw0ZzttsubhddFDN/XHklfPe7UZn/5jebriS/+ipceGH83I8/PnrN1123+fcygwMOiF7vf/4zrkBssEHMvFHYm+0eV+H//e9I8N94IwZs/vjHbfteG1HFt3HFV60OIiJNU8U3PbNmwYEHwpgxcZl4333bnvTm5ef0TfJS86RJUckt9fdh9dUjWSpnhodFi+Cgg+IS+hVXtP54bLZZJL+/+AWcfDJ88knx55XT5pBXiYrvggWRPH796xF/UlZZJZa3vv/+aB9ZY43onx40KKrK+asEzz4bz9tll0h0X3sNLr645aS3UIcOMRjzpZfgq1+NVf4OPDD+3Xrr+HvffPNY8e+uu6ICveGGyX2vOar4Nu7xVcVXRKRpkydHP18l1ELiu2hRVDCL3d59N+aXPemkllcTe/rpSCyOOALOPTdWxkvS8OGR7D31VMSShPxsDuXMhHT66VF53WefSICb61H95JPoT11nnahINp6urVxbbBHJ3ciRsOOOUX3cbLMVn1NOm0PhfmfNghkzVly5rtDUqdEOMGBA+XG7x4wNgwZFy0Za1l8/9v+Tn0R7yBVXxN/6JpvE7/Jpp8VViGKDCcvRuXP8TRxzTLRJ9OgRCe4GG0DPnol8K82p78RXFV8RkdJ99llcftxkk8q8XxZbHZYujRa5O++M27//HVWv9dZbftthh6hirbVW9DMee2wkPSedFKPqC/9zd4/Baj/7WVRkDzkknbg7dIjBcKNGJZP4ukfie8MN5b1unXXiw9Pll8el7333jdaOxnMAf/JJVB7XWgv+9re2J715PXvGvLdXXgm77gq//nW0KOST9/Hj45J/OTp0iGWhx49fsQXhzTfjvW64IT4I7bwz3HNP+TG/8kpcabn77spMt9qxY3ww2WefKA4+/XTMC1y4OEoSunWL/vMKq+/Et0+fGJG4bFn84qriKyLStFdfhf79K7fIT48eUVHNn6OrZcECeOCBWJ3snntg7bXj8vtVV0WPYnOx7bBDVDkffTQS27PPjtkITjopKoX//d8x08H48S3PjNBWxx8fq6j97ndtr6xNnRofhHbYofzXdu4M3/52fAi48MLoNf3Wt6K62LNnXOI++OC4KnvNNcklvXlmcOKJMGRIVNgffjgu63/+eUz/1ZpVCYcNiz7foUNjtoUbbogPRV/9akyntuWWUVFdtCgSvnLcc0+0OSR9HErR0BDtCO1Ifff4du4ciW5+OUdVfEVEmlbJ/l6IhLJbt5jGqhoeeSSqkuusExXKwYNjcYPnn4dzzonL5aUk5B06xOX90aPh5ZcjwT3yyFhowSwqyGknvRBXOffZJwa5tdUtt0QC35YK5GqrxXF84YWoLG68cSw4ccgh0Lt3OklvoS22iGpmjx6RwP/lL/Ezbc175lfP23RTeOaZWEb43Xcj6R02DHr1iuWbWzOjwT33tLvks5rKTnzNrK+ZjSv4uqeZ3WtmD5jZbWbWJdkQU1bY56uKr4hI0yZPbl01rC2q0ef70ktR0T3xRDjqqEhg7rsvZkRYf/227btvXzjjjJhT9ZlnIrkrtwLYFkkNcmtptbZy5Ht4H3ooKusNDXDddcn3ORfTrVscj3POicFaQ4e2bj9DhsQVgffeg2uvjQVGOnde8TkHHhjtCuVYsCD6svfcs3VxyUrKSnzNrBdwDVC4jt9RwMXuvg8wA9gvufAqoLDPVxVfEckIM+tkZm+Z2ZjcbSszu9LMnjSzswuet9K21FS64guV7fP94IMYMb/rrjEY7KWXYgBOGgWRDh2iV7oSPZuF9tgjKujPPNP6fbz2WhSMdt45ubggPlTdfXckjpVIegsdcUS08px2Wuten6/qd+3a9HPyiW85HzoeeiiS6rYOKJP/KLfiuxQ4AvjPWcjdL3H3B3NfNgAfFnuhmY00s4lmNnFm4UwK1aaKr4hk09bAje4+3N2HAwOBju4+BBhgZgPN7NDG21KNqBqJbyUqvosWxUwKm28e/csvvxx9uUkP5smCDh3avpLbrbdGD241+67T0LdvugnmZptFG8WUKaW/Jt/fK4lp9rfWzC4rqDaMAU5193lNPHcI0Mvdiy4I7u6j3H2Quw9qaGhoc+CJUcVXRLJpJ2CEmT1tZlcCewGjc489AAwFhhfZlo6FC2O6phTm1WxW2hXfm26Kyutzz8Ul5d/9LuYybc+OPx5uvrn1vdM33hhTr0l5zMprd3BX4puCZhNfdz8pX23I3X5Z7Hlm1hv4E3BCGkGmShVfEcmmZ4C93H0w0BnYH3g399gcoC/RdtZ420oSueI2deryilUlpVnxfe21aG34xz8iAd5oo3TeJ2v69o0BV7fcUv5rX3ghlrrdbbfk46oH5SS+zz8fPcgbb5xuTHWmzdcpcoPZbgLOdPfpbQ+pwlTxFZFsesHd38/dnwj0AfINhD2I8/fCIttWksgVt2q0OUC6Fd+f/AS+973yFyxoD447LgbWleu662IasvbW5lApw4dHQjtnTsvPVbU3FUn85n4T2B44K9cScUQC+6wcVXxFJJuuM7NtzKwjcDBwCstbGbYBpgGTimxLR7US37QqvpMmwdixsTxqPTrooEjAppdRr1q6NOanPeaY9OJq77p2jWr5/fe3/Ny771bim4JWDZvMDbTI378UuDSpgCpOFV8RyaZfAjcABtwB3A6MM7O1ibaHnQAvsi0dU6bEHLCVttpq6SS+Z54ZFd/u3Vt+bnu0yirRp/v3v8NZZ5X2mkceicU7Gi/zK+XJtzt87WtNP2f27Jg+UC0lidO1inzF9/PPYcmSyq1IJCLSDHef4u5bu/tW7n6Wu88nBrNNAHZ393nFtqUW0JQplZ/DF6Lim3Srw8MPx9LLJ56Y7H5rzXHHxdRhpU6vde21qvYm4YADYl7opUubfs4DD0RbhHKSxCnxzVd8FyyIykKl51QUESmRu89199HuPqO5bal45plYaKDSkm51cI9q7znnrLzAQL3Zccf496mnWn7uwoVw552x4py0Tb9+UTlv7rirzSE1SnzXWCOWLP7oI/X3iog0Zd11q1MYSHpw2y23xNU9TccVP89jjy1tkNutt8ZMEGuumX5c9aC52R2WLo2KsBLfVCjx7dQJevaMBn/194qIZEuSFd8lS6Kf9fzzNStB3jHHwOjR8NlnzT/vuuvU5pCk5hLfZ56BL34xKsOSOP3lQ/T5/vvfqviKiGRNkhXfv/0t2jWqMUgvq/r1g222iTaGprzzTsyCcdBBlYurvdtpJ3j77Ti2janNIVVKfCEu3fz736r4iohkTVIV30WL4Be/gAsu0FiOxvKD3Jpyww1w2GExFZcko1Mn2HffmKu3sXvuiYqwpEKJL6jiKyKSVUlVfP/85xjMNXhw2/fV3hx2WMxpnJ/as5B7JMXHHlv5uNq7Aw9cOfF9//2YcWTIkOrEVAeU+IIqviIiWZVExXfuXLjoIjj33GRiam969IAvfzkqu4099xx8/HF9rm6Xtv32g0cfXbG/+t57Ye+9NeNIipT4QlR8X39dFV8RkaxZbbVY3vWf/4Snn45510uddzbv17+Ggw+GTTdNJ8b2oKl2h/ygNg0GTN4aa8RqiI89tnyb2hxS16qV29qdNddcPo+viIhkxxe+AD/8Idx0E7z5Ztw+/xz694cNNojbeutFEtH41rs3fPABXH55LM8rTdt9d5g1K1YLyy9UsmRJVIHHjatubO1ZfnaHffaJ3+uHHoK//KXaUbVrSnwhKr6giq+ISNaYxaC0QvPmwbRpyxPhd96BqVNjmdfC20cfxSCi006LeYilaR06wNFHR9X3ooti24MPxgeMgQOrGlq7duCB0WP9+9/D+PFxrPv2rXZU7ZoSX1g+IbcqviIi2dezZ0zBtc02zT9v2bJIkldfvTJx1brjjovK7/nnxwcGDWpL39ZbR4/vq6+qzaFC1LQDqviKiLRHHTpAr16avqxUm2wS8/o++GDMpHHvvXDEEdWOqn0zizl7775b8/dWiCq+oIqviIgIRIX32mtjWq3dd49eaUnXgQdGH/ucOTBoULWjaffKrviaWV8zW6nTPbf92WTCqrDeveNTlyq+IiJSz448Miq9l16qJYorZY89old9//01e0YFlHWEzawXcA3QvcjDvwFqc1mXjh2hTx9VfEWkZpnZlWb2pJmdXe1YpIb17g177hlz26vftDJ69Ihjfeih1Y6kLpTb6rAUOAL4V+FGM9sD+BiY0dQLzWwkMBKgX79+Zb5tBRx/fIxeFRGpMWZ2KNDR3YeY2VVmNtDdX6t2XFKjfvjDaHNYZZVqR1I/br212hHUjWYTXzO7DNikYNMj7v5LKxgoYGZdgJ8AhwC3N7Uvdx8FjAIYNGhQmbOPV8CFF1Y7AhGR1hoOjM7dfwAYCvwn8c184UGyZccd4ybSDjWb+Lr7SSXs4wzgEnf/yDRyVkSkGroD7+buzwG2L3ww84UHEZEKSaKLei/gFDMbA2xrZlcksE8RESndQpaPseiBpqoUESmqzdOZufuu+ftmNsbdT2zrPkVEpCyTiPaGCcA2wCvVDUdEJJtalfi6+/BytouISKpuB8aZ2drA/sBOVY5HRCSTdDlMRKTGuft8YoDbBGB3d59X3YhERLJJK7eJiLQD7j6X5TM7iIhIEeZe+QG+ZjYTmN6Kl/YBZiUcTlIUW/myGhcottbIalxQvdjWd/eGKrxvk9pw/oXs/oyzGhdkN7asxgWKrTWyGhdk7PxblcS3tcxsortnciFrxVa+rMYFiq01shoXZDu2WpLV45jVuCC7sWU1LlBsrZHVuCB7sanHV0RERETqghJfEREREakLtZb4jqp2AM1QbOXLalyg2Fojq3FBtmOrJVk9jlmNC7IbW1bjAsXWGlmNCzIWW031+IqIiIiItFatVXxFRERERFpFia+ItJqZ9Tazvc2sT7VjKZTVuEREkpLV81xW48qrmcTXzK40syfN7Oxqx1LIzDqZ2VtmNiZ326raMQGYWV8zG1fwdSaOX2FcWTl2ZtbTzO41swfM7DYz65Kh41Ustqofs1xsvYC7gMHAo2bWkIXj1kRcmThmtSoLP9disnIOaUzn37Ji0vm3dbHp/NtKNZH4mtmhQEd3HwIMMLOB1Y6pwNbAje4+PHebXO2Acr941wDdc19n4vg1jovsHLujgIvdfR9gBnAkGTheTcR2Btk4ZhA/v9Pc/VzgfmAPsnHcGsd1Atk5ZjUnK+ePJmTlHPIfOv+WTeff1tH5t5VqIvEl1qDPL8X5ADC0eqGsZCdghJk9nfu0lYVloJcCRwDzc18PJxvHr3FcmTh27n6Juz+Y+7IBOJpsHK9isS0hA8csF9tj7j7BzHYlPt3vSwaOW5G4PiEjx6xGDScDP9cmZOIc0ojOv2XQ+bfVsen820q1kvh2B97N3Z8D9K1iLI09A+zl7oOBzsABVY4Hd5/v7vMKNmXi+BWJK1PHzsyGAL2At8nA8SpUENuDZOuYGfGf6VzAychxaxTXs2TomNWgTJw/mpCpcwjo/NtaOv+2Ki6df1uhVhLfhUDX3P0eZCvuF9z9/dz9iUCWLgPmZfX4ZebYmVlv4E/EZZlMHa9GsWXmmAF4OAV4AdiZjBy3RnGtnaVjVoMy9ffQSKb+HpqQ1eOXmWOn82/r6PzbOln5A2zJJJaX7bcBplUvlJVcZ2bbmFlH4GDg+WoHVERWj18mjp2ZdQFuAs509+lk6HgViS0TxywX24/M7Njcl6sDF5CB41Ykrr9m5ZjVqMz8PRSRmb+HZmT1+GXi2On82+rYdP5tpZpYwMLMVgPGAQ8D+wM7NbpkUzVmtiVwA2DAHe5+VpVD+g8zG+Puw7N2/AriysSxM7OTgfNY/gf5N+A0MnC8isT2KHAYGfh9yw2WGQ2sAkwBzgTGUuXjViSuS4HrycAxq0VZO38Uyso5pBidf0uOR+ff1sWm828r1UTiC/85mHsDY919RrXjqTU6fuXR8WodHbf2ST/XttHxK4+OV+vouJWmZhJfEREREZG2qJUeXxERERGRNlHiKyIiIiJ1QYmviIiIiNQFJb4iIiIiUheU+IqIiIhIXVDiKyIiIiJ1QYmviIiIiNQFJb4iIiIiUheU+Io0YmZ9qh2DiIiIJE+Jr1SVmXVp9HVDiu+1pZl9rYXndAYeNbNm1xM3s9XM7Awzs9zXXcysQ8HjnQq/FhHJmtacf83sAjM7pcz32dzMfmpmHYs89h0zW9XMHjGzrczsB7nz62Vmtms57yNSCv3HLFWTOwk+ZWZbF2y+28z2LnjOLmZ2jJkdbWbHmtkeZrbIzCaa2Xtm9mb+fu75DWb2j4LbjWbW0cw6Ae8CB+We1ymftDbyF+BOYA8zO7RIzNua2ePAPcC5wITc1xcAj5vZLDN7ARgPbJfEcRIRSVop59/c8/qb2fUFmz4FPit4fJSZrZe7b2b2We6c/KmZbZ572g+Ak4HhRULpBJwNLAG6A0e6+3xgD+CdJmJ/yMyeMrMxjW5PmdkDZRwGqUNKfKVq3H0pcCbwXQAz2wgwd3+wyNMvzP27GHjV3QcBVwG/AHbKbQdYBegH/Ag4DdgR2B6YCDwEbGxmzwOTgE3yOzezDmb2B6CDu/8YOAT4HzP7eaOqSA/gdSJB3tnddwReA37o7jsDbwK7uvuO7j6p9UdHRCQ9ZZ5/129mV5sDHXP7dOBtYBfgBXd/0cyGARsCg4CLzWzD/AvNbFXgPWAm0ItIjMfmKr093P2N3PNWafSen9G0xc08JkKnagcg9cvMxhCf8D/LVU3J3X8V+Je7/8DdxwPjzewsd78213+7Ue75/Yjqw4mA516/FPgcOAH4EPjc3Z8xsx8B67r7lWZ2FfCku7+ci6M/cDkwH/iuma2V29dxwGXAv3OvuYxIoGcDjwP7ALj7N8zsx2Z2GrA68LqZTXD3EYkfNBGRBJRy/m3lrpcAuwH35M6tlwEHuPu7ZvYd4CEz+7673wJ0AXYFBhOJ8UvAHKI63MHMJhLn+QVmtq27L8i9x0Qiwf200XuvCnRuZdxSJ5T4SjV9Dhzv7tMKN5rZ8cDA3P0vECezDmbWl7hK8bq7DzWzc4jqUaBdjQAAIABJREFU699z/8LyBPhzclUBMxsAfAJ8ycyeAvYFLs2dlHsSbQtnAYcBV+Te+2PiMltH4ADgf4DViMrEJUSV4rdmth1RHe4I3AI8CXTL7UtEJKtaPP/mrAJsa2bP5b5eC/g8l8QCbMTKucSXgWeAx4CXgTvNbClxnrwCuMjMlhDn3jeAdYgk9iNgLrAz8Gd3/5WZ/QP4rbsvMLODgdOBeayc9OZ1NrNxwKm66ibFKPGVaupCVAU+b7S9N/DP3P3jgJ/ktk0GVuq7zenQ6N8lgOVumwGbAq8AewG/IaoML7r7vWa2kbt/AlwNkE+o3f3qgv1/K/fYrsDzuW19cs97xMx2BsYBhwO/aiZOEZEsKOX8CzAd6JlrjcDMfg684+5XNLPvR4CDgSPcfUKuVaGbu881s37ARe7+uZltAMwCfgj8EvhZ7nWdgR1y++oH/BvA3W8HbjezK4grb4uBvrnnfZB7XW93/0p5h0LqSVUT31wF72Z3H9bE4ycDR+S+XB14yt1PqlR8ki53H56/b2ZbufvkIs/5M/BnM3vZ3Tc1szWBJblLYGsTVd3vAC/md5X7tydxIuzg7neb2YFA4SCO8e7+29x7fGJmk3KvXZbfb66i0RHo6u6b5p47NvcfxU+AZ4EhZnYZ8H5uv/0L7ouIZFIp59+cq4kq7UPFHsxViLd39/8p2PetZnYBcJ6Z/Yo4V+6aazl7hThPfkAk2d8jxmp0IAYWLyYqvmPMrBfwBXef0+httycqvp8T1WKIwctdcjeRJlVtcFvuF/oaoseoKHe/1N2H5/5AxxF9mNIOmNneZranmX3LzLYCLsidQDGzDc3scis+HdhMYKfCwW25+wfmHs+f9NYkTqyzcrM3bEaMFh4KfJ8VL+UB7EdcTmu8378D/+nVzV1qux3YABhKJMZjWV5dNmAhGjgqIhlV5vl3Se7WlHwCWrj/brntqxBtDxOIZHYb4D53/wAg14owlGhd2JgYxPxK7grcrcS59pEi79kp99g/iIHKk3L3b0fnXmlBNX9BlhLV3PkQfyhmdrOZjTWzvxQ+0czWAfq6+8QqxCnpOIpIUncEGohWgsNzj71B9MleWvD8rmZ2BjFbw6RcxfcE4Ge5+5Ny/cDTiZ7bYcTJ8ZjcSOOljd5/WaOvOxE9u+vlN5hZd2Ianv+MIM5davsucEPuA9m03G0sUXkYD+yOPqSJSHaVe/4t13nA08SVuGlEAjyOmLXhX42euzYxWO0+YiBcvi/3VqIl7c7CJ+cKGafn9v0O0RM8N3f/ReCMNsQtdaBqrQ65efqw5VOpjgSmuPvPzexWM9va3V/IPXYKbfsjlAzJJZS7AyeRm1fX3d8jBpFBnPBOIS6TfZ24FDaeGCRxn7tfkNtPsV5czOxo4qT7IXCHme1OVGLvNrPFxFWGFwtf4+7v51oWji/YfAJwrbu/3ehbcOBbZrYfUTkeRQx2Oxg4kuUzQfzT3Rsn2CIiVVPm+Xd1YnDxFWa2MPd448FtqwM3F7xFJ3c/NZegrkIMDH7P3U83s0eBr+emMfssV5R4najW/pIYZLy1xRzqZxNX50aZ2XHuPi63/wHAOcSAZVje6rBZ7t/OZvaUu89rw2GSdixLg9s2AXY2s+HEH9I6wAu5yy27E38Q0j50Ak5z98/M7DXgKjPLT1TeDZiROxEfb7GS0APuPgvAYlGKgUTyWdiLa8B1RB/aN4E93X22mZ3J8iluDnT3d8xsJwp+n3IJ79Dcl4uJ37/FxJy9c81s/1xcG+cGeHQBLnf3c8zsaiLh3Rv4mrtPM7OvECfyXsQADBGRrCj5/Jvb9tXmdmZmRxJtDPlqbH51tt7ErA1LAC8oct1EnJO/bGYfAXcDzwGHuvuMXOHi20Rr2qu58Rc/zbWZdSd6jmexfAafhY3+NeAuMzvM3T8s/bBIvbD4wFXFAMzGuPtwMzsVmOfufzOzEcCb7j7VzHYDDnH3U6saqGSCxWISi73RL67FKkQd3H2xmXVy9+Z60pKOyYBVc31pIiIiklFZagK/HNjfzMYC/02s/gIx5+rYqkUlmeLunzdOenPbl7r74tz9iiW9ufdzJb0iIiLZV/WKr4iIiIhIJWSpx1dEpG4Vm7ecOEdvDtzt7ufknndlKdtERGRlVUl8+/Tp4/3796/GW4uIVMykSZNmuXtDKc9190vJzV5jZn8i2r02d/chZnZVblDnVkDHlra5+2tNvY/OvyJSD5o6/1Yl8e3fvz8TJ2pKXhFp38xseitesw6xDKsDo3ObHyBmHtmuxG1NJr46/4pIPWjq/JulwW0iIrJ83vLuxDKsAHOIZLjUbSsws5FmNtHMJs6cOTPF0EVEsk2Jr4hIRhTMWz6GmJe0a+6hHsT5utRtK3D3Ue4+yN0HNTSU1HkhItIuKfEVEcmOYcBTuSn7JrF8YZVtiKVfS90mIiJFaFYHEZHsKJy3/HZgnJmtDewP7ET0/ZayTUREilDFV0QkI9z9x+5+a+7+fGA4MAHY3d3nlbqtGrGLiNSC9pP4vv8+3HZbtaMQEUmMu89199HuPqPcbTXllltg1qxqRyEidaD9JL6PPw4nnwxLl1Y7EhERKcf558Mzz1Q7ChGpA+0n8Z07Fz74ACZMqHYkIiJSjpkzYcGCakchInWgfSW+XbrArbdWOxIRESnHrFkwf361oxCROtC+Et9DDok+X/dqRyMiIqVYtChuqviKSAW0r8R3+PC4//zzVQ1FRERKlB/UpsRXRCqgfSW+vXrBoYdqdgcRkVqRX0JZrQ4iUgHtM/FVn6+ISG3IJ76q+IpIBbS/xHenneLS2WuvVTsiERFpyaxZ0LGjKr4iUhHtL/Ht0AEOPljtDiIitWDmTOjXTxVfEamI9pf4QszuoHYHEZHsmzULBgxQxVdEKqJ9JL7LlsVJc/XV4+vhw+HVV+Hdd6saloiItGDmzEh8VfEVkQpoH4nv/PnQvXv0iUEsZDFiBNx+e3XjEhGR5qniKyIVlHjia2aXmNlBSe+3WYVtDnn5xSxERCS7VPEVkQpKNPE1s2HAWu5+Z5L7bVGxxHfffeGZZ2D27IqGIiIiZchXfJX4ikgFJJb4mlln4HJgmpl9Jan9lqRY4tutG+y5J9x1V0VDERFpi8KrZmZ2pZk9aWZnFzxe0raakZ/VYfHiuImIpCjJiu+xwIvAhcBgM/tu4YNmNtLMJprZxJn5CcuTUizxBS1mISI1pfCqmZkdCnR09yHAADMbWOq2Kn4L5Vm6FD76CNZYA77wBVV9RSR1SSa+2wGj3H0G8Hdg98IH3X2Uuw9y90ENDQ0Jvi1NJ74HHgiPPgoLFyb7fiIiCSty1Ww4MPr/27v38Kiqqw3g7zKECAEJCCLRAiLY1oooRgUVDChatBZLtVrrXUtFq7b1SrW1tfbmvfppK4KKiBb8rPgJVVEB5apcRG71LiBUMMgdIQRY3x8rY5JhMjkzs88t8/6eZx4mZybn7JwMJ+usvfbe1S9PAnBCBtviYd06m42noMACXw5wIyKfuQx8PwLQpfp5GYDlDvedXn2Bb+vWQO/ewMsvB9YUIqIs1ek1A3AVgMScjOsAtAdQ7HFbHb72uOWiogJo29ae77MPM75E5DuXge9IAP1E5E0AVwK42+G+06sv8AVY7kBEcZHca/YmgGbVr7WAXa+3eNxWh689brlYuxZItIcZXyIKgLPAV1U3q+rZqtpXVXuranCrR6QLfAcNAl56CaisDKw5RERZSO4164yasoUeAJYBmOdxWzww40tEAWsSdgOcSBf47r8/8J3vAJMnAwMHBtsuIiLvRgJ4TETOBVAIq939PxEpBTAQQC8ACmCah23xkJzxZeBLRD5rHCu3pQt8AS5mQUSRl6LXbDks+J0NoJ+qblTVTV62hfMTZCE548tSByLyWf4Evi+8YFPnEBHFhKquV9Vx1XW/GW2LhYoKZnyJKFD5Efh26QKUlgIzZgTXJiIiSm/tWmZ8iShQ+RH4AsAppwBvvhlMe4iIqGHM+BJRwOIf+O7eDWzcaJOgp/PNbwIffhhMm4iIqGGczoyIAhb/wHfzZqB5c6BJAxNUdOsGfPBBMG0iIqKGcTozIgpY/ANfL2UOAHDIIcz4EhFFhSozvkQUuPwJfPfbD9ixw9aGJyKicG3dCohYjx3AjC8RBSJ/Al8RZn2JiKKidrYX4OA2IgpE/gS+AOt8iYiionZ9L8DpzIgoEPkV+DLjS0QUDbWnMgOY8SWiQORX4MuMLxFRNNRevAJgxpeIApFfgS8zvkRE0ZCc8W3WDKiqsgcRkU/yK/BNZHxV/W0TERGll5zxFQFatGC5AxH5Kr8C39atgb33Btas8bdNRESUXnLGF+CUZkTkOyeBr4g0EZEVIjK1+tHdxX49ySTwBVjnS0QUBcnTmQEc4EZEvnOV8T0cwDOqWl79WORovw3LNPBlnS8RRUyq5IGIjBSRWSJya633edoWC8nTmQEc4EZEvnMV+PYC8D0Rebv6Itwk+Q0iMkRE5orI3IqKCkeHRXYZXwa+RBQtdZIHALoBKFDV3gC6iEg3ERnsZVtoP0GmmPElohC4CnznADhZVY8BUAjgtOQ3qOpwVS1T1bJ2yRe7XGST8WWpAxFFS53kAYCTAYyrfm0SgBMAlHvcFg/M+BJRCFwFvgtV9fPq53Nh2Qr/qQIbNgAlJd6/hxlfIoqe5OTBQACrql9bB6A9gGKP2/bgW49btqqqLMBNTlow40tEPnMV+I4WkR4iUgDgTADvOtpveps329yPhYXev6drV+Cjj4Ddu/1rFxFRZpKTB20BNKv+ugXsWr3F47Y9+Nbjlq116yzoLSiou71lS2Z8ichXrgLf2wGMBrAAwCxVfc3RftPLtMwBsHki27QBVq70p01ERJlLTh5chZqyhR4AlgGY53Fb9KWaygzgdGZE5Ls9BqFlQ1UXwwZnBCubwBeomdKsY0f3bSIiytztAJ4GIAD+D8B4ANNEpBRW9tALgHrcFn3Ji1cktGxprxER+STeC1hkG/hySjMiihBVXayqh6tqd1W9RVU3wQauzQbQT1U3et0Wzk+QIWZ8iSgkTjK+ock140tEFFGquh41MzZktC3yUk1lBnBwGxH5jhlfIiIKVqqpzABOZ0ZEvsvPwJcZXyKi8DDjS0Qhyc/A9+CDgRUrbC5JIiIKFjO+RBSS/Ax8i4qADh2AZcucN4mIiBpQ3+A2ZnyJyGf5GfgCrPMlIgpLuunMmPElIh/lb+DLOl8ionBwOjMiCkn+Br7M+BIRBU+1/oxv8+ZAZSWwc2fw7SKivBD/wLekJLvv7daNgS8RUdA2bwYKC4FmzfZ8TYR1vkTkq/gHvrlkfFnqQEQUrPqyvQkMfInIR/ENfFWBDRuyD3w7dQJWrwa2b3fbLiIiql999b0JnNKMiHwU38B3yxagaVN7ZKNJE6BzZ+Djj502i4iI0qhv8YoEZnyJyEfxDXxzKXNIYJ0vEVGw6lu8IoEZXyLyUX4HvqzzJSIKVkOlDsz4EpGPnAa+ItJeRN5xuc96MeNLRBQ/Xga3MeNLRD5xnfG9G0CKOWp8wIwvEVH8eBncxowvEfnEWeArIv0BbAWwup7Xh4jIXBGZW1FRkfsBmfElokamdq+ZiIwUkVkicmut1z1tizROZ0ZEIXIS+IpIUwC/AXBzfe9R1eGqWqaqZe3S3e175SLwPeAAmxKNF1kiioa7ATQTkcEAClS1N4AuItLN67YQ2+4NpzMjohC5yvjeDOBhVd3gaH8NcxH47rUX0LUr8NFHbtpERJSlpF6zcgDjql+aBOCEDLal2rfbHrdcMONLRCFyFfieDOAqEZkK4AgRGeFov/VzEfgCVu7AOl8iClGKXrNiAKuqn68D0D6DbXtw3uOWC2Z8iShETVzsRFX7Jp6LyFRVvdzFftNyFfgecgjrfIkobF/3mokIAGxBzUDhFrAkhddt0bVjhy0+VFJS/3uY8SUiHzm/SKpquet9psSMLxE1HnV6zQCcgZqyhR4AlgGY53FbdH35JbDvvlZmVh9mfInIR04yvqFwmfF99NHc90NElKXkXjMA3wcwTURKAQwE0AuAetwWXQ2VOQDM+BKRr6LdLZaOy4wvSx2IKCJUtVxVN8EGrs0G0E9VN3rdFk6rPWpoYBvABSyIyFfM+O63H1BVBaxbB7Rpk/v+iIgcUNX1qJmxIaNtkeUl48sFLIjIR/HM+Kq6C3xFmPUlIgqC14wvA18i8kk8A9+tW4HCQqCoyM3+uHQxEZH/vGR8i4uB7duBnTuDaRMR5ZV4Br6usr0JzPgSEfnPS8ZXBGjRwqY9IyJyjIEvwIwvEVEQvGR8AU5pRkS+YeALMONLRBSEioqGM74A63yJyDcMfIGaRSxU3e2TiIjqWruWGV8iChUDX8D2tffewJo17vZJRER1eS11YMaXiHzCwDehc2dgxQq3+yQiIqNqSxaz1IGIQsTAN6G0FPjvf93uk4iIzMaN1rPmZRpKljoQkU8Y+CYw8CUi8o+XqcwSmPElIp8w8E3o0IGBLxGRX7zW9wLM+BKRb5wGviLSRkQGiIjH2/osMeNLRBQvzPgSUQQ4C3xFpDWACQCOATBFRDze2mfBr8D388/d7pOIiAwzvkQUAU0c7utwAL9S1dnVQXBPAK843H8NZnyJqBESkTYAjgLwjqquDbs9TnldvAJgxpeIfOMs46uqb1QHvX1hWd9Zrva9Bwa+RNTIpOo1E5GRIjJLRG6t9T5P2yLH6+IVgAW+zPhS1Jx6KrB8editoBy5rvEVAOcAWA+gKum1ISIyV0TmVlRUZH8QVX8C37Ztbbqdykq3+yUi8ibRa/ZHWG9ZfwAFqtobQBcR6SYig71sC+0nSCeTjO8++zDjS9GyZg0waZI9KNacBr5qrgKwEMD3k14brqplqlrWzutdfypffQUUFNh8kC7ttRew//7A6tVu90tE5EGKXrNTAYyrfnkSgBMAlHvcVoezxEMuMs34MvClKJkxAygsBKZODbsllCOXg9tuEpELq78sAbDB1b7r8CPbm8ByByIKUVKvmQJYVf3SOgDtARR73FaHs8RDLji4jeJs+nTgooss8FUNuzWUA5cZ3+EALhCRNwEUwDIP7jHwJaJGKqnX7DgAzapfagG7Xm/xuC16OJ0Zxdn06cAFF1iP84cfht0ayoHLwW3rVXWAqvZV1StVfbol8jPw5SIWRBSSFL1mf0FN2UIPAMsAzPO4LXqY8aW42roVWLoUOOYYoLyc5Q4x53I6s2Aw40tEjdNwAONE5HIAiwGMB/CmiJQCGAigF6z8YZqHbdFSWQls2wa0auXt/cXFwPbtwK5dlmEjCtPbbwOHH25ji/r1A159FRgyJOxWUZai2SWWjt+BLxexIKIQpOg12wgbuDYbQD9V3aiqm7xsC+cnSCNR5iDi7f0iQIsWLHegaJg+HTihulOlvByYMoV1vjHGwLc2ZnyJKEKqg+Fxqro6022RkslUZgms86WoqB34du4MNG0KfPBBqE2i7DHwrY2BLxGRe5lMZZbAwJeiYNcuYPZs4Ljj7GsRK3eYMiXcdlHWGPjWxsCXiMi9bDK+HOBGUbBokcUGtT+/HOAWaxzcVlubNrZAxrZtQLNmDb+fiIiAjz8GBg60QWklJXs+Fi0C2u8xvXB6zPhSFNQuc0goLwduvtnqfL3WrVNkMPCtTcSmNPv8c6BLF3+OQUTU2MyaBRxyCHD77cCGDXs+OnQAfvjDzPbJjC9FwfTpdlNXW+fOlhx77z3g298OpVmUPQa+yRJz+TLwJSLyZuFCq4Hs2dPdPpnxpbCpWuD7xz/u+Vqi3IGBb+ywxjcZ63yJiDLz7rs2z6lLzPhS2FassMFtqRJhrPONLQa+yTiXLxFRZhYudB/4MuNLYZs+HTj++NR1vInAl/P5xk68Al9VZnyJiKLkiy9sQPA3vuF2v8z4UthSDWxL6NTJBnP+5z/BtolyFq/Ad9s2u/Pyc8YFBr5ERN4tWmTZXtej25nxja/33gN27w67FblLF/gCNau4UazEK/D1O9sLMPAlIsrEwoVAjx7u98vAN562bbPPw4MPht2S3KxfDyxfDhxxRP3v6dePdb4xxMA3GQNfIiLv/KjvBVjqEFezZlkZwB132KDHuJo1Czj6aKBJmsmvTjwReOMN1vnGDAPfZAx8iYi88yvwZcY3nqZOBc46C7j3XuDHP7ZFoeKooTIHAOjY0T6nS5YE0yZywlngKyKtROQlEZkkIs+LSFNX+/5aEIFvq1ZAVRWwZYu/xyEiirudO21wz3e+437fzPjG09SpVvt6/vnAkUcC110Xdouy4yXwBVjuEEMuM74/AXCvqp4CYDWA7zrctwki8K29ehsREdXvgw+AAw4AWrRwv29mfOPnq6+A+fNtMRMR4OGHgVdeAZ5/PuyWZaay0n6OXr0afi/n840dZ4Gvqj6sqq9Wf9kOwBeu9v21IAJfgHP5ElHgUvWaichIEZklIrfWep+nbYHwq8wBYMY3jmbPts9D4kaoVStgzBjgiiuAlSvDbVsm5s0DvvlNu/lqSCLwbQyzWOQJ5zW+ItIbQGtVnZ20fYiIzBWRuRUVFdntPMjAl3W+RBSs5F6zcwEUqGpvAF1EpJuIDPayLbAW+xn4MuMbP1OnWtd/bb17A1dfDVx4oa2CFgczZtjCFV4ceKDFJazzjQ2nga+ItAHwIIBLk19T1eGqWqaqZe3atcvuAAx8iaiRStFrdj6AcdVfTwJwAoByj9vqcJJ4SMWvqcwAWxxg27b4BEtkc9qWl++5fdgwqwe/667Am5QVr/W9CZzPN1ZcDm5rCuBZAMNUdbmr/dbBwJeIGrlErxmAzwCsqt68DkB7AMUet9XhJPGQip8Z3732suCXA43j4auvgHfesfreZAUFwFNPAffdB7z9dvBty8Tu3ZllfAHW+caMy4zvZQB6ArhFRKaKyDkO920Y+BJRI5bUa7YFQGKZyhaw67XXbf5bv94enTv7dwyWO8THrFmW/S8uTv16x47AQw8B550X7d/p++9bffkBB3j/nvJym8+Xdb6x4HJw299VtbWqllc/xrra99cY+BJRI5Wi12weasoWegBYlsE2/y1aBHTvbplZv3CAW3wkpjFL56yzrAb45z8PokXZmTEjszIHwILkffe1/xMUeWmWJImg9euBkhL/j8PAl4iCV7vX7BYAjwO4QERKAQwE0AuAApjmYZv/3n3XvzKHBGZ842PqVOC3v234ffffb8sAP/888IMf+N6sjE2fnlmZQ0Ki3MGvmndyhiu3pZIIfLkMIREFJEWv2SjYwLXZAPqp6kZV3eRlWyAN9rO+N4EZ3+x98UVwf8PS1fcmKy4GnngCuPJKwOVAS8A+k+edB0yblv0+Mh3YltCvHwe4xQQD31QSc/cx00BEIVLV9ao6TlVXZ7rNd0EEvvmc8Z00yTKi991nq+N5CWI//RS4806grAxo3x544QX/2wk0XN+b7PjjgQsuAIYOdRecjxkDnHQScNBBtlTyRRcBa9Zkto/PPwfWrQO+/e3Mj3/yycBbb9nMFUyaRVp8At+qKusead48mONxEQsiotR27bJ5S7t39/c4+Zrxfe45W/L3lFOA994DTj3VArorrrBgtvbNQCLYPfpo4NhjgU8+Af76V2D0aOBvfwumvV7qe5PdfrsF9P/8Z27H3rEDuOYa4LbbgNdfB/74R9tvu3bAYYfZ6nFep8RLzOaQTd16u3YW+I4dC5x7LmcjibD4BL6FhbYqjEgwx2OdLxFRap98Yn/oW7Xy9zj5mPEdNcoGf73yimVEH3kEWL4c+Pe/gW7dgP/5H/v71L9/3WD3L3+xv1n/+IdlPs85x5aUXrjQ/zZPmbLnwhUN2Xtv+1mvvTb7v7X//a+dh08/BebMqemBaNkSuPtua9fYsXaO6ptGbfVqu9H41a9svuFMf47aOna0MovmzW3hjo8+yn5f5Jt4DW4LEgNfIqLUgihzACzjm0+B70MPWbZ2yhTgW9+q2S4CHHqoPa67zrKJU6da8FheDjRJ8ae8sNDqaB94ABgxwr82b90KLFhggV6mysosuP/pT4EJEzJLbE2bZpnVn/0MuPXW1Fnaww6z8zRmDHDmmcAZZwCXXw7Mn2/Z3RkzrISyd2/L9I4Ykd3AttqaNQMee8xuQI4/Hnj8ceC003LbJzkVn4xv0Bj4EhGlFlTg27JlfEodtm+3AV6jRllwOmCATXM1aBDw5psN133++c/AvffafLC1g95UWrQAvvc9qytNFfQmDBli2cy1azP/ebyaNcvKEL3W9ya75RYrK3z8cW/vV7USjrPOskD1t79NX5ogYmUjS5cCRUVW+ztrFtCnD/Dii3ZuJk4Efv1r4MQT059Pr0QsoP/Xvyyo/8Mf/J3j9667gLZtrbb52WdZZtEABr71YeBLRJRakIFvVDO+u3ZZfeqPfwx85zs28PrCC21QWrt2wC9/aVnJgQMty3jssdbtvnNn3f2oWhf7U0/Z+w86yF0b27WzAXKPPupun8myqe+trWlTu1m46SYr6Uhn/nz7eZ54woLXgQO9H6ekxLLfS5fa9//0p5ZB93Me6uOPtxKMl18GBg8GNvow4coDD1h2+dVXrUxj5EiLXwYNsp/zyy/dHzPmGPjWh4EvEVFqQczhC0RzcNuOHRZcfOtbVm976qnA008DGzbYAgZjxgA332zd21262IC0996zzOZDDwFdu9pctps3Wxbw6qstaHnjDfu749o119gAr6oq9/sGcg98ARsked11wKWX7pkZVQVeesnqlgcNsmBy5kw7t3FQWmqlK/vsA1x/vdt9Dx9uvQSTJwNHHmkZ/pdfBlasAM4+2zLaXbrYuRsxIro3kQFj4FsfBr5ERHvatMmmiera1f9jRSnju20b8OCD9nOPHWuBxLRpwMUX21ReRUX1f+9ee9WUPDz7rA3UPuggCxjffdeKan5iAAAXdklEQVRmI2jb1p92H3GEBT/PP+9+37nU9ya7/nqbD/jvf7evKyut/KF7d8uIX3IJ8PHHwA03WB1tnDRtajdDkye72+eoUVZC8dprQKdOdV8rKbHyjueeszKSn//cyjk6drTeh9mz3U+5pmrHmjHDei9uv93+b5SX28DBlSvdHi8HHNxWHwa+RER7WrzYuvYLCvw/VhQyvps2WTB2//1WrvDcczabQraOPtpKJJYts5kaLr7Y/2k6r7nGMoM/+pHb/eZa31tbkyYWzB13nM208NhjFvTef79lLIOa0ckv3/62fZZWrgQOPDC3fY0dazcDkyc3fAPavLmVh/zgB3ZeR42yOZSLiiwIPv/87G+6li2zm8GXXrLnLVvaDV2XLvZvnz7AT35iWejDD7cBhjfckN08yQ4x41ufDh3s7oUTURMR1QiqvhcIP+P70kvAwQdbVnbSJGD8+NyC3to6d7ZZF4KYm37QIAu45s51u18XZQ61HXKIDdRatcrO/csv2wC+uAe9gP0MffrktqocYJn7a6+16e4aGgSZbP/9rZb6gw+s7GbePAuczznHgul16xreh6qVmpx9NnDUUfZzPf20rRS4Zo1lk59+2uZTvuwyG+R5zz02tdtBB9kAwsGDbc7jdKqqLMM/f35mP6MHDHzrU1xs08H4UYxORBRXQQa+YWZ8R4607vUXXrA/5H4v1uGnJk2Aq66ygVAuZTN/b0MuucSyvUF9xoLUt6+Vu2Tr3/+26dsmTszt8yhiAejo0TYHcnm5Pe/c2cpWfv97C2BrL/xRVWU9Fb16Wca4Tx/L8t59t2X9W7RIf8w2bYDf/MaO16+fBdv9+tnN5PjxFhwPHWqLthx8sO3vpJMso+wYSx3SSZQ7lJSE3RIiomhYuNB9l3l9wsj4qtof/ieftCDlkEOCPb5fLr/cAorVqy3zl6utWy0T7qK+N1/07Zv9DBuvv25lMS+8YJlWV1q3toBz6FCrq54+3TLtP/2pxT8DBtjCKaNGWcZ22DCbDznbUqfiYhvQecUVlmW+7z67we3a1Uqovv99+5x27my10T5g4JtOIvA99NCwW0JEFD5Vm7kgqOxn0IFvVZUFAAsWWP1q+/bBHdtvbdrYDcvw4Tb3ba5mzrSZBIIo1WgsevSwMo61azOrq/3qKystGD/e3xuNoiLLsp50kpWcrFxpJT5Lltixe/Z0d6zCQqsvPv98d/v0yGmpg4i0F5EcC1gihAPciIhqJAaw7LtvMMdr0cL+6Ps5+X/Cli1WC/vf/1rtamMKehOuvtoG6u3Ykfu+XNf35oOCAgtcp0/P7PsmT7bSj759/WlXfQ480KaYu+cet0FvyJwFviLSGsAoAA6Gd0YEA18iohpB1vcCNg1Y8+b+r0S1Zo0FcR06WFdyQ/WKcXXYYdad/Oyzue+LgW92sqnzffFFW6mPnHCZ8d0F4BwAEZttPAcMfIkoQMm9ZiIyUkRmicitmW7zxcKF1l0bJL8HuH3wgU2hdcYZNjdvYaF/x4qCa66xJX9zmbGI9b3ZyzTwVQUmTLDPJznhLPBV1U2qWu8UCCIyRETmisjciooKV4f1FwNfIgpIcq+ZiAwGUKCqvQF0EZFuXrf51sigM76Av3W+O3bYdFnDhgG33dY4ps1qyOmn2zK2DU0nlQ7re7NXVmYr+Xm9mVuwwM5zYxlkGQGBTWemqsNVtUxVy9q1axfUYXPToQMDXyIKSnKvWTmAcdXPJwE4IYNtdThLPIQR+PqZ8X3mGZsL9fLL/dl/FBUU2Epef/tb9vuYMoVlDtkqKrK5oGfO9Pb+F1+0bG8+3JQFhPP4plNaaotYEBH5LEWvWTGAVdXP1wFon8G25H3nnnj46ivgs8+Czzz5lfFVtTlIr7/e/b6j7tJLbYGItWsz/15V4F//Ak47zX278kUm5Q4TJrC+1zEGvulw9TYiCs8WAM2qn7eAXa+9bnNvyRLgm98MvgZ2n338CXwnTbIs2oAB7vcdda1aWcnDM89k/r1z5tjCBr16uW9XvvAa+K5eDXz4IXDCHp04lAPnF0hVLXe9z9A0a2a1NV6W8SMicmseasoWegBYlsE29959N5zVtFq29KfUIZHtzdcu5IsuskUJMvXkk8CFF+bveXOhVy+r3d22Lf37Jk60lcx8WsghX3EBi4YkBrgFNW8lUX127wb+8Q+bW/HEEy1rQ43ZeADTRKQUwEAAvQCox23uhVHfC/hT6rBgAfCf/wDnnut2v3Fy0knWo7lkiU1x5sWOHbba1ttv+9u2xq642BaBeeut9LXSEyYAgwcH1qx8wVKHhnBmB4qKu+4CHnnE1i4/8EDg2GNtNPprr1n9JTUKiV4zVd0EG7g2G0A/Vd3odZsvDQtjKjPAn8Ft99xj03rlcyatoMBWzXrySe/f8+9/20qmBx3kX7vyRUPlDtu328IVAwcG16Y8wYxvQxj4UhS89Zb9sZ4zB+jUyS6Ks2fbhfG226wb+uijbV7NAw4A9t/fVp7af397NNYJ+Rs5VV2PmhkbMtrm3EcfNY6M72efWRfygw+622dcXXSR1Tj/6U8WCDdk9Gjgggv8b1c+6NsXuO+++l+fOtWywpksbUyeMOPbEAa+8VNVBTzwAPDpp2G3xI2NG4HzzrOlRjt1sm17721dZLffDsyYYV2WN95o25csAcaMAW66CTj1VGC//Szw7doVePzxUH8UirHly+2zFDTXGd8HHgAuvhgoKXG3z7g69FD7G/faaw2/d906e9/ZZ/vfrnxw/PGW0KiqSv06Z3PwDTO+DSkttVowio9hw4BXXrGg8KSTbADL0UeH3arsqAJXXGFZmR/+sP73tWxpXWKpusVUbcnXBQusXuzMM4HWrf1rMzVOXjKCfnCZ8d24EXjsMWD+fDf7awwuusjKHU49Nf37xo616wvHFrhRUgIcfLB9Fo89tu5ridXaJk4Mp22NHDO+DUlMaUbxMH488L//a91En35qXf9nnWWDwV580QaIxckTTwCLFqXvEmuIiAUPffoAP/gB8Oc/O2seke9cZnxHjLAAL9FzQjbAb+LEhs9xYjYHcqe+Ot/Fi+26feihwbcpDzDwbQhLHeLjk0+AIUMsM7Hvvhbs/eIXwMcfA0OHAr/7nY1eHjHCamT94mre5/ffB264AfjnP21qPRd+9ztg5EircySKA1cZ36oq4P77geuuy31fjUnbtkC/fsCzz9b/ng8+sETCKacE1658UF/gmyhz4JRxvmDg2xAGvvGwfbtldm+9dc9uoyZNLKsxd67VyT7/vI1KvvNOd12oqraM5/e/b3M/f+97Vk/75ZfZ7a+y0tp8xx3AYYe5aSNgn+crrgB++1t3+2zMVDmPd9hcLWAxbhzQrRtw1FG576uxaWhO39GjgR//2K6l5E6fPsD06bYgSG2JZYrJFwx8G7L//rZ6Sty6yMOkCixdahfS1auDOeYvfmGDt66+uv73iNiAsIkTrQZ4wQKgSxebFSHbAHXHDusC7NkTuPJKWw3p009tMNqECbb/k0+2gDuTkpmbbrLg/Gc/y65d6dx4o01LtGiR+303FhUVwL332k3H0KFhtya/uVjAIp+XJ/bitNNsLMsnn+z52u7dwFNPsczBD+3b22Px4pptFRU2QPnEE8NrVyPH27eGFBVZMf/ateGMaN62DXjnHeum27Ur9aO01AKvTJYS3b3bAp+FC20Z0u7dc+tO/+wz4PXXax5NmwJHHAH88peWKbjxRv/q6saMsWm95s713jV0+OHA00/bcpB33mmZoEsusW7Q0tKGv3/tWptT96GHrHziT3+y2sG9qu8lzzvPHl99Bbz8MvDcc8Cvf23vPe004JhjgLKy1CPLJ0ywrPQ77/jT1dWqlQ0AHDbMjuW33bvtxmL3bsuGN28e3kCpdHbtsmVsR4600euDBtkNS58+Ybcsv7kodZg82W5Sv/tdN21qbJo2tR6m0aMtEVDb9Ok2K8wRR4TTtsYuUe6QmCP7pZdsUHZRUbjtasREXdUjZqCsrEznzp0b+HGz1qOHZS+D/I8/f77Voo4da5m/RLBQ+9GkiQVan3wCLFtmA7lOPNGymmVldSdnV7Wa0cmT7TF1qtXBHnGE1W+9/74Ffz171jx69LALXmJWgPXrax4bNlgX8Lx5FiR8+SXQv7/9hz35ZMt0igBr1lhd3fDh1nVz883At77l7jwtXWo/8+uv5zbH6MqVNk/uqFHAj35kwc62bRa4bttW9/natRbMDh5smebu3b0do7Ky5sZgzhz7HR9wgAXBxxxjM0/stx9w3HFWb+dnwFVZab+HJ57IPbOwebNNy7NkiWX416yxfxPPv/jCgpcmTYCtW+0cNm1qqxcVF9tnO/G8vkdhoX0GN22y423aVPNIBEX77pv60aaNfY6Limy6t6Kius+3brWboCeesOzLZZfZzZqD0esiMk9Vy3LekUOxu/5u2mS/l0cfBTp3thvo0tLMbp4GDrRpuC691Ldmxt7cucA559h8zbVvuC+/3JIjN9wQXtsas6eeAl54oabG+uyzLTlyySXhtqsRqO/6y8DXi4EDrQv9tNP8Pc769Za9HDnSnl96qc032bFjw9/75ZfAtGnAG2/Y48MPrdb1uOMsMJ482YKH/v3t0a+frf6VsH27dbfMn1/zWLzYssAbN1qQ0Lq1ZShbt655fthhFuj26FGT7UxlwwbLjj7wgAV0v/61Bde52LLFAsbrr3f3B62iwia2/+gjC8iaNbNH7ectWlhJQ/v2uR1r504L3OfMsSVA58yxLPzvfgfccouTHyetMWPsZ501y3tmWdU+TzNn2vfNnGmftcSNUocONYtmJBbQ2G+/PW/CEjcSW7fW/JvusWOHBc/77FPzqP317t12I/bll3s+1q2zY2zfbgF/8r8iNtvFZZc5X6CBga8DqlbrvnSpzSW8fLndfJaW1gTCBxxQ90anTZua56tW2bV72TJm0dJRtev5I48AJ5xg27Zts3O7eLG3njDK3IoVlvRYvdp6dvfbzxJRuf59IQa+ObnsMsumXn65+30nBkWNGGF1l9/9rh2nf//0gWRDNmywLqpZs+yPQ//+NVlYr6qqbD+tWrlb2nPrVsvc3H23BUsPP1w3APdK1VYQKixsXIsy7NwZ3ACS3butZ+CWW9LPEQxY0HHbbdYl16SJ3VAlHkcemd9Lv6bBwNcnlZVWXrV8uQW0q1ZZsiBxo5P4d906237PPcC114bd6ui78067kX30Ufv6n/+06+srr4Tbrsauc2c7xytXWlLorbfCblGjUN/1lzW+Xvgxs0NlJfDMMzaAZtcuG8T04IOWpXChpMRmFshl5ZfCQqBdOzftSSgutvKAoUPtItuzJ/C3v1nXsleqNq/twoW2bG9jEuSo6b32Av76V+Cqq2w2ilQ14jt32u/pvvtsxox77gG+8Q1Os0PhKiqywaxduzb8XlV+Xr06/3zL+j7wgPVuPfkklygOQqLOd8kSzuYQAM7q4IXLRSzWrrVuu86d7W767rutG+maa9wFvXFQVAT85jdWyP+HP9jAioamjVK1rHhZmdVFPfeclSBQ9gYMsK7ikSP3fG3hQiuXmTrV6v+uvdbKbhhEUJzw8+pdaal1u48fb13vM2daGRD5q08fC3xffJHLFAeAga8XLjK+771n86d262Zdc6++agOkTjklvy/MRx1lA+Q6dLD6ypdfTv2+yZNtbfMbb7Su+Xnz7FxS7v76V1veecsW+7qqqma55yuvtC44rnRFlB8SSxg/84wtb15cHHaLGr++fe1mY8eOmtkdyDdO+1VFZCSAQwFMVNU7XO47VKWlNthpyhQbQZ4YVZ54vnmz1a5WVtoHt7Ky7mPrVguchw61AJhF63U1a2Zd6WecYSNZTz8duOsuu+DOnGmZ4RUrgN//3kYdR3EqrDjr2dNmAkn8Di6+2D7z77yTXf01EcXXmWcCP/+5zfaTqPUlfx1ySM3CR/mcCAuIs8BXRAYDKFDV3iLymIh0U9UPXe0/VN26Wc3s739vI8kTj8TI8k6dLEgrKrJBPonpkmo/cp0nNx/07w+8+66VfRx5pNXvLV1qq4xdeCFXDfLTHXdYAPzgg3bTceGFvADHTKNNPFCwmje3qRonTbIbYvKfiM1MdPrpYbckL7iMJMoBjKt+PgnACQC+DnxFZAiAIQDQ0cv0XFHSujUwY0bYrcgPJSXWzfbCCzYH7PPPcwqiIHTpYjXn3bvb9EUUK4068UDBGzbMantzmVmIMvPnP4fdgrzhMvAtBrCq+vk6AHUmaVXV4QCGAzadjsPjUmM0aFDYLcg/XNUqzsrRWBMPFLyDD7YHUSPk8nZuC4BEX34Lx/smIqL6JSce6gwkUNXhqlqmqmXtXE9RSEQUIy6D03mwLAMA9ACwzOG+iYiofkw8EBF54LLUYTyAaSJSCmAggF4O901ERPVLJB5mwxIP74fbHCKiaHIW+KrqJhEpBzAAwJ2qutHVvomIKC0mHoiIPHDaHaaq61V1nKqudrlfIiKqn6pugg1wmw2gHxMPRESpcWJUIqJGQFXXo2ZmByIiSkFUg59ZTEQqACzP4lvbAljruDmusG2Zi2q7ALYtG1FtFxBe2zqpaqSmUcjh+gtE93cc1XYB0W1bVNsFsG3ZiGq7gIhdf0MJfLMlInNVtSzsdqTCtmUuqu0C2LZsRLVdQLTbFidRPY9RbRcQ3bZFtV0A25aNqLYLiF7bOOUNEREREeUFBr5ERERElBfiFvgOD7sBabBtmYtquwC2LRtRbRcQ7bbFSVTPY1TbBUS3bVFtF8C2ZSOq7QIi1rZY1fgSEREREWUrbhlfIiIiIqKsMPAloqyJSBsRGSAibcNuS21RbRcRkStRvc5FtV0JsQl8RWSkiMwSkVvDbkttItJERFaIyNTqR/ew2wQAItJeRKbV+joS5692u6Jy7kSklYi8JCKTROR5EWkaofOVqm2hn7PqtrUGMAHAMQCmiEi7KJy3etoViXMWV1H4vaYSlWtIMl5/M2oTr7/ZtY3X3yzFIvAVkcEAClS1N4AuItIt7DbVcjiAZ1S1vPqxKOwGVX/wRgEorv46EucvuV2Izrn7CYB7VfUUAKsBnIsInK962nYzonHOAPv9/UpV/wjgFQD9EY3zltyuSxGdcxY7Ubl+1CMq15Cv8fqbMV5/s8Prb5ZiEfjC1qBPLMU5CcAJ4TVlD70AfE9E3q6+24rCMtC7AJwDYFP11+WIxvlLblckzp2qPqyqr1Z/2Q7A+YjG+UrVtp2IwDmrbtsbqjpbRPrC7u5PRQTOW4p2bUNEzllMlSMCv9d6ROIakoTX3wzw+pt123j9zVJcAt9iAKuqn68D0D7EtiSbA+BkVT0GQCGA00JuD1R1k6purLUpEucvRbside5EpDeA1gA+QwTOV2212vYqonXOBPbHdD0ARUTOW1K73kGEzlkMReL6UY9IXUMAXn+zxetvVu3i9TcLcQl8twBoVv28BaLV7oWq+nn187kAotQNmBDV8xeZcycibQA8COuWidT5SmpbZM4ZAKi5CsBCAMchIuctqV2lUTpnMRSp/w9JIvX/oR5RPX+ROXe8/maH19/sROU/YEPmoSZt3wPAsvCasofRItJDRAoAnAng3bAblEJUz18kzp2INAXwLIBhqrocETpfKdoWiXNW3babROTC6i9LAPwFEThvKdr1j6ics5iKzP+HFCLz/yGNqJ6/SJw7Xn+zbhuvv1mKxQIWIrIPgGkAXgcwEECvpC6b0IjIYQCeBiAA/k9Vbwm5SV8TkamqWh6181erXZE4dyIyFMCfUPMf8nEAv0IEzleKtk0B8ENE4PNWPVhmHIAiAIsBDAPwJkI+byna9XcAYxCBcxZHUbt+1BaVa0gqvP56bg+vv9m1jdffLMUi8AW+PpkDALypqqvDbk/c8PxlhucrOzxvjRN/r7nh+csMz1d2eN68iU3gS0RERESUi7jU+BIRERER5YSBLxERERHlBQa+RERERJQXGPgSERERUV5g4EtEREREeeH/ASuTYzTcxjz7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 画出4个时间序列的图像\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2,figsize=(10,6))\n",
    "for i, ax in enumerate(axes.flatten()):\n",
    "    data = df[df.columns[i]]\n",
    "    ax.plot(data, color='red', linewidth=1)\n",
    "    ax.set_title(df.columns[i])\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ccee775a-f4da-4495-9709-291affef9898",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>V1_x</th>\n",
       "      <th>V2_x</th>\n",
       "      <th>总理赔金额_x</th>\n",
       "      <th>总出险数量_x</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>V1_y</th>\n",
       "      <td>1.0000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>V2_y</th>\n",
       "      <td>0.0000</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>总理赔金额_y</th>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>1.0000</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>总出险数量_y</th>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>1.0000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          V1_x   V2_x  总理赔金额_x  总出险数量_x\n",
       "V1_y    1.0000 0.0000   0.0000   0.0000\n",
       "V2_y    0.0000 1.0000   0.0000   0.0000\n",
       "总理赔金额_y 0.0000 0.0000   1.0000   0.0000\n",
       "总出险数量_y 0.0000 0.0000   0.0000   1.0000"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# step2：Granger’s Causality Test ， 检验不同序列之间存在互相影响\n",
    "maxlag=12\n",
    "test='ssr_chi2test'\n",
    "variables=df.columns\n",
    "def grangers_causation_matrix(data, variables, test='ssr_chi2test', verbose=False):    \n",
    "    df = pd.DataFrame(np.zeros((len(variables), len(variables))), columns=variables, index=variables)\n",
    "    for c in df.columns:\n",
    "        for r in df.index:\n",
    "            test_result = grangercausalitytests(data[[r, c]], maxlag=maxlag, verbose=False)\n",
    "            p_values = [round(test_result[i+1][0][test][1],4) for i in range(maxlag)]\n",
    "            min_p_value = np.min(p_values)\n",
    "            df.loc[r, c] = min_p_value\n",
    "    df.columns = [var + '_x' for var in variables]\n",
    "    df.index = [var + '_y' for var in variables]\n",
    "    return df\n",
    "\n",
    "grangers_causation_matrix(df, variables = df.columns)   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d3af795-4feb-4224-ada2-4c8bb03dafeb",
   "metadata": {},
   "source": [
    "在输出结果中，index以y结尾，表示响应变量，column以x结尾，表示预测变量，如果p值小于0.05表明存在格兰杰因果性。  \n",
    "因此根据检验数据，完全有理由使用VAR模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "69d3d619-dec0-427a-a8c0-d7329cf4a1cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# step3：ADF测试，检验单个变量是否平稳\n",
    "def adfuller_test(series, signif=0.05, name='', verbose=False):\n",
    "    r = adfuller(series, autolag='AIC')\n",
    "    output = {'test_statistic':round(r[0], 4), 'pvalue':round(r[1], 4), 'n_lags':round(r[2], 4), 'n_obs':r[3]}\n",
    "    p_value = output['pvalue'] \n",
    "    def adjust(val, length= 6): return str(val).ljust(length)\n",
    "\n",
    "    # Print Summary\n",
    "    print(f'    Augmented Dickey-Fuller Test on \"{name}\"', \"\\n   \", '-'*47)\n",
    "    print(f' Null Hypothesis: Data has unit root. Non-Stationary.')\n",
    "    print(f' Significance Level    = {signif}')\n",
    "    print(f' Test Statistic        = {output[\"test_statistic\"]}')\n",
    "    print(f' No. Lags Chosen       = {output[\"n_lags\"]}')\n",
    "\n",
    "    for key,val in r[4].items():\n",
    "        print(f' Critical value {adjust(key)} = {round(val, 3)}')\n",
    "\n",
    "    if p_value <= signif:\n",
    "        print(f\" => P-Value = {p_value}. Rejecting Null Hypothesis.\")\n",
    "        print(f\" => Series is Stationary.\")\n",
    "    else:\n",
    "        print(f\" => P-Value = {p_value}. Weak evidence to reject the Null Hypothesis.\")\n",
    "        print(f\" => Series is Non-Stationary.\")    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "02817b4a-1b37-4875-b4d1-2908043693e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Augmented Dickey-Fuller Test on \"V1\" \n",
      "    -----------------------------------------------\n",
      " Null Hypothesis: Data has unit root. Non-Stationary.\n",
      " Significance Level    = 0.05\n",
      " Test Statistic        = -4.4212\n",
      " No. Lags Chosen       = 2\n",
      " Critical value 1%     = -3.633\n",
      " Critical value 5%     = -2.949\n",
      " Critical value 10%    = -2.613\n",
      " => P-Value = 0.0003. Rejecting Null Hypothesis.\n",
      " => Series is Stationary.\n",
      "\n",
      "\n",
      "    Augmented Dickey-Fuller Test on \"V2\" \n",
      "    -----------------------------------------------\n",
      " Null Hypothesis: Data has unit root. Non-Stationary.\n",
      " Significance Level    = 0.05\n",
      " Test Statistic        = -4.4355\n",
      " No. Lags Chosen       = 2\n",
      " Critical value 1%     = -3.633\n",
      " Critical value 5%     = -2.949\n",
      " Critical value 10%    = -2.613\n",
      " => P-Value = 0.0003. Rejecting Null Hypothesis.\n",
      " => Series is Stationary.\n",
      "\n",
      "\n",
      "    Augmented Dickey-Fuller Test on \"总理赔金额\" \n",
      "    -----------------------------------------------\n",
      " Null Hypothesis: Data has unit root. Non-Stationary.\n",
      " Significance Level    = 0.05\n",
      " Test Statistic        = -11.3894\n",
      " No. Lags Chosen       = 1\n",
      " Critical value 1%     = -3.627\n",
      " Critical value 5%     = -2.946\n",
      " Critical value 10%    = -2.612\n",
      " => P-Value = 0.0. Rejecting Null Hypothesis.\n",
      " => Series is Stationary.\n",
      "\n",
      "\n",
      "    Augmented Dickey-Fuller Test on \"总出险数量\" \n",
      "    -----------------------------------------------\n",
      " Null Hypothesis: Data has unit root. Non-Stationary.\n",
      " Significance Level    = 0.05\n",
      " Test Statistic        = -4.0456\n",
      " No. Lags Chosen       = 1\n",
      " Critical value 1%     = -3.627\n",
      " Critical value 5%     = -2.946\n",
      " Critical value 10%    = -2.612\n",
      " => P-Value = 0.0012. Rejecting Null Hypothesis.\n",
      " => Series is Stationary.\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for name, column in df.iteritems():\n",
    "    adfuller_test(column, name=column.name)\n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c8d1d6-3da7-4ab3-9b5d-a068a0ed6ba2",
   "metadata": {},
   "source": [
    "#### 无需差分，均已平稳"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "dd261ec2-8d57-4a55-9826-7bd531500511",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Name   ::  Test Stat > C(95%)    =>   Signif  \n",
      " ----------------------------------------\n",
      "V1     ::  71.73     > 40.1749   =>   True\n",
      "V2     ::  18.92     > 24.2761   =>   False\n",
      "总理赔金额  ::  6.96      > 12.3212   =>   False\n",
      "总出险数量  ::  0.0       > 4.1296    =>   False\n"
     ]
    }
   ],
   "source": [
    "# step4: 协整检验，检验多变量平稳性\n",
    "def cointegration_test(df, alpha=0.05): \n",
    "    out = coint_johansen(df,-1,5)\n",
    "    d = {'0.90':0, '0.95':1, '0.99':2}\n",
    "    traces = out.lr1\n",
    "    cvts = out.cvt[:, d[str(1-alpha)]]\n",
    "    def adjust(val, length= 6): return str(val).ljust(length)\n",
    "\n",
    "    # Summary\n",
    "    print('Name   ::  Test Stat > C(95%)    =>   Signif  \\n', '--'*20)\n",
    "    for col, trace, cvt in zip(df.columns, traces, cvts):\n",
    "        print(adjust(col), ':: ', adjust(round(trace,2), 9), \">\", adjust(cvt, 8), ' =>  ' , trace > cvt)\n",
    "\n",
    "cointegration_test(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ed7b96d6-0323-4851-912a-13f741aa3fdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# step5：划分训练集和测试集\n",
    "nobs = 2  # 最后四个时间点作为测试集\n",
    "df_train, df_test = df[0:-nobs], df[-nobs:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "00c28215-4ff2-4a86-bbab-d5c942a2793c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# # step6：使用VAR之前，先差分处理使单个变量变得平稳\n",
    "# df_differenced = df_train.diff().dropna()\n",
    "# for name, column in df_differenced.iteritems():\n",
    "#     adfuller_test(column, name=column.name)\n",
    "#     print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "420af00b-dab3-4346-bb1c-f6b2ef83decb",
   "metadata": {},
   "source": [
    "#### 所有序列已经平稳，无需进一步差分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "8e0c5eec-c659-4507-b8a5-00aaf5164c54",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lag Order = 1\n",
      "AIC :  46.90844303468482 \n",
      "\n",
      "Lag Order = 2\n",
      "AIC :  42.13550046778399 \n",
      "\n",
      "Lag Order = 3\n",
      "AIC :  41.51288762727691 \n",
      "\n",
      "Lag Order = 4\n",
      "AIC :  40.391389179791155 \n",
      "\n",
      "Lag Order = 5\n",
      "AIC :  39.20407281729118 \n",
      "\n",
      "Lag Order = 6\n",
      "AIC :  36.40375409911314 \n",
      "\n",
      "Lag Order = 7\n",
      "AIC :  -18.791001865185827 \n",
      "\n",
      "Lag Order = 8\n",
      "AIC :  -83.46012502567673 \n",
      "\n",
      "Lag Order = 9\n",
      "AIC :  -85.58173313913319 \n",
      "\n",
      "Lag Order = 10\n",
      "AIC :  -148.6390578872253 \n",
      "\n",
      "Lag Order = 11\n",
      "AIC :  -143.42246805918305 \n",
      "\n",
      "Lag Order = 12\n",
      "AIC :  -140.83983541818804 \n",
      "\n",
      "Lag Order = 13\n",
      "AIC :  -145.9659805198062 \n",
      "\n",
      "Lag Order = 14\n",
      "AIC :  -145.19420367389625 \n",
      "\n",
      "Lag Order = 15\n",
      "AIC :  -147.32843146544747 \n",
      "\n",
      "Lag Order = 16\n",
      "AIC :  -153.78727960015098 \n",
      "\n",
      "Lag Order = 17\n",
      "AIC :  -155.4313815298087 \n",
      "\n",
      "Lag Order = 18\n",
      "AIC :  -161.7643913088828 \n",
      "\n",
      "Lag Order = 19\n",
      "AIC :  -179.36158718675955 \n",
      "\n"
     ]
    }
   ],
   "source": [
    "# step7：选择模型阶数并训练，根据AIC值，lag=4时达到局部最优\n",
    "# df_differenced = df_differenced + 0.0001\n",
    "model = VAR(df)\n",
    "for i in range(1,20):\n",
    "    result = model.fit(i)\n",
    "    print('Lag Order =', i)\n",
    "    print('AIC : ', result.aic, '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "f08f8d87-1915-4ce6-bf3c-452198eb4cc5",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  Summary of Regression Results   \n",
       "==================================\n",
       "Model:                         VAR\n",
       "Method:                        OLS\n",
       "Date:           Sun, 06, Nov, 2022\n",
       "Time:                     23:02:09\n",
       "--------------------------------------------------------------------\n",
       "No. of Equations:         4.00000    BIC:                    43.4441\n",
       "Nobs:                     34.0000    HQIC:                   41.4325\n",
       "Log likelihood:          -811.629    FPE:                5.16493e+17\n",
       "AIC:                      40.3914    Det(Omega_mle):     1.02023e+17\n",
       "--------------------------------------------------------------------\n",
       "Results for equation V1\n",
       "===========================================================================\n",
       "              coefficient       std. error           t-stat            prob\n",
       "---------------------------------------------------------------------------\n",
       "const           -2.563204         0.748769           -3.423           0.001\n",
       "L1.V1            9.386370         8.322012            1.128           0.259\n",
       "L1.V2            0.036654         0.031050            1.180           0.238\n",
       "L1.总理赔金额         0.000001         0.000000            2.478           0.013\n",
       "L1.总出险数量        -0.002204         0.001302           -1.693           0.091\n",
       "L2.V1          -25.901622        12.192844           -2.124           0.034\n",
       "L2.V2           -0.097874         0.044639           -2.193           0.028\n",
       "L2.总理赔金额         0.000000         0.000000            0.321           0.748\n",
       "L2.总出险数量         0.000216         0.001837            0.118           0.906\n",
       "L3.V1           10.944574        11.029488            0.992           0.321\n",
       "L3.V2            0.043671         0.040259            1.085           0.278\n",
       "L3.总理赔金额        -0.000000         0.000000           -1.576           0.115\n",
       "L3.总出险数量        -0.000389         0.001781           -0.218           0.827\n",
       "L4.V1           -3.150727         7.398943           -0.426           0.670\n",
       "L4.V2           -0.016372         0.027694           -0.591           0.554\n",
       "L4.总理赔金额         0.000000         0.000000            1.793           0.073\n",
       "L4.总出险数量         0.001558         0.001106            1.409           0.159\n",
       "===========================================================================\n",
       "\n",
       "Results for equation V2\n",
       "===========================================================================\n",
       "              coefficient       std. error           t-stat            prob\n",
       "---------------------------------------------------------------------------\n",
       "const          673.862349       198.405632            3.396           0.001\n",
       "L1.V1        -2235.578726      2205.131945           -1.014           0.311\n",
       "L1.V2           -8.761856         8.227631           -1.065           0.287\n",
       "L1.总理赔金额        -0.000152         0.000060           -2.515           0.012\n",
       "L1.总出险数量         0.590633         0.344943            1.712           0.087\n",
       "L2.V1         6781.294050      3230.808936            2.099           0.036\n",
       "L2.V2           25.619445        11.828378            2.166           0.030\n",
       "L2.总理赔金额        -0.000016         0.000062           -0.264           0.792\n",
       "L2.总出险数量        -0.063428         0.486790           -0.130           0.896\n",
       "L3.V1        -2917.084321      2922.547606           -0.998           0.318\n",
       "L3.V2          -11.604716        10.667592           -1.088           0.277\n",
       "L3.总理赔金额         0.000077         0.000050            1.532           0.125\n",
       "L3.总出险数量         0.094800         0.471887            0.201           0.841\n",
       "L4.V1          883.849037      1960.541065            0.451           0.652\n",
       "L4.V2            4.504061         7.338295            0.614           0.539\n",
       "L4.总理赔金额        -0.000045         0.000025           -1.766           0.077\n",
       "L4.总出险数量        -0.404456         0.292993           -1.380           0.167\n",
       "===========================================================================\n",
       "\n",
       "Results for equation 总理赔金额\n",
       "===========================================================================\n",
       "              coefficient       std. error           t-stat            prob\n",
       "---------------------------------------------------------------------------\n",
       "const      2651954.544438    861559.588385            3.078           0.002\n",
       "L1.V1      6369187.343754   9575597.969843            0.665           0.506\n",
       "L1.V2        16901.252316     35727.789148            0.473           0.636\n",
       "L1.总理赔金额        -0.211080         0.262162           -0.805           0.421\n",
       "L1.总出险数量      5076.303177      1497.886893            3.389           0.001\n",
       "L2.V1     20829939.330880  14029513.091231            1.485           0.138\n",
       "L2.V2        83319.488878     51363.723633            1.622           0.105\n",
       "L2.总理赔金额        -0.514681         0.268184           -1.919           0.055\n",
       "L2.总出险数量     -1817.300475      2113.845790           -0.860           0.390\n",
       "L3.V1     -3912751.049411  12690914.478754           -0.308           0.758\n",
       "L3.V2       -22994.769003     46323.111108           -0.496           0.620\n",
       "L3.总理赔金额         0.607556         0.218536            2.780           0.005\n",
       "L3.总出险数量      2765.577472      2049.129053            1.350           0.177\n",
       "L4.V1     -1515793.377415   8513482.870331           -0.178           0.859\n",
       "L4.V2         4232.121612     31865.920246            0.133           0.894\n",
       "L4.总理赔金额        -0.286042         0.110221           -2.595           0.009\n",
       "L4.总出险数量     -4043.001221      1272.296315           -3.178           0.001\n",
       "===========================================================================\n",
       "\n",
       "Results for equation 总出险数量\n",
       "===========================================================================\n",
       "              coefficient       std. error           t-stat            prob\n",
       "---------------------------------------------------------------------------\n",
       "const         1328.381816       333.401398            3.984           0.000\n",
       "L1.V1        -2490.449052      3705.510088           -0.672           0.502\n",
       "L1.V2          -11.891328        13.825735           -0.860           0.390\n",
       "L1.总理赔金额        -0.000291         0.000101           -2.873           0.004\n",
       "L1.总出险数量         2.196347         0.579644            3.789           0.000\n",
       "L2.V1        11883.093995      5429.060665            2.189           0.029\n",
       "L2.V2           45.689734        19.876440            2.299           0.022\n",
       "L2.总理赔金额        -0.000064         0.000104           -0.615           0.539\n",
       "L2.总出险数量        -0.538280         0.818004           -0.658           0.511\n",
       "L3.V1        -3860.452110      4911.057437           -0.786           0.432\n",
       "L3.V2          -16.064318        17.925852           -0.896           0.370\n",
       "L3.总理赔金额         0.000183         0.000085            2.167           0.030\n",
       "L3.总出险数量         0.281836         0.792960            0.355           0.722\n",
       "L4.V1          119.168804      3294.498867            0.036           0.971\n",
       "L4.V2            2.810154        12.331291            0.228           0.820\n",
       "L4.总理赔金额        -0.000100         0.000043           -2.356           0.018\n",
       "L4.总出险数量        -0.814532         0.492346           -1.654           0.098\n",
       "===========================================================================\n",
       "\n",
       "Correlation matrix of residuals\n",
       "               V1        V2     总理赔金额     总出险数量\n",
       "V1       1.000000 -0.999580 -0.708120 -0.927830\n",
       "V2      -0.999580  1.000000  0.707685  0.926951\n",
       "总理赔金额   -0.708120  0.707685  1.000000  0.628070\n",
       "总出险数量   -0.927830  0.926951  0.628070  1.000000\n",
       "\n"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 选择lag=4拟合模型\n",
    "model_fitted = model.fit(4)\n",
    "model_fitted.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "e55619f2-df25-4cbe-9092-1536138fe18b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "V1     : 2.75\n",
      "V2     : 2.75\n",
      "总理赔金额  : 2.53\n",
      "总出险数量  : 2.76\n"
     ]
    }
   ],
   "source": [
    "# step8：durbin watson test，检验残差项中是否还存在相关性，这一步的目的是确保模型已经解释了数据中所有的方差和模式\n",
    "out = durbin_watson(model_fitted.resid)\n",
    "def adjust(val, length= 6): return str(val).ljust(length)\n",
    "for col, val in zip(df.columns, out):\n",
    "    print(adjust(col), ':', round(val, 2))  # 检验值越接近2，说明模型越好"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "1a94d6ea-6395-482b-bc84-42e8ad680820",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>V1_2d</th>\n",
       "      <th>V2_2d</th>\n",
       "      <th>总理赔金额_2d</th>\n",
       "      <th>总出险数量_2d</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>-3.1301</td>\n",
       "      <td>824.0439</td>\n",
       "      <td>3052220.3584</td>\n",
       "      <td>1559.8748</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>-3.1559</td>\n",
       "      <td>828.7594</td>\n",
       "      <td>3736476.9421</td>\n",
       "      <td>1764.4167</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     V1_2d    V2_2d     总理赔金额_2d  总出险数量_2d\n",
       "36 -3.1301 824.0439 3052220.3584 1559.8748\n",
       "37 -3.1559 828.7594 3736476.9421 1764.4167"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# step9：模型已经足够使用了，下一步进行预测\n",
    "lag_order = model_fitted.k_ar\n",
    "forecast_input = df_differenced.values[-lag_order:]\n",
    "fc = model_fitted.forecast(y=forecast_input, steps=nobs)\n",
    "df_forecast = pd.DataFrame(fc, index=df.index[-nobs:], columns=df.columns + '_2d')\n",
    "df_forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "08dc14be-e6a1-4d7f-9a7a-77811e2298b6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>V1</th>\n",
       "      <th>V2</th>\n",
       "      <th>总理赔金额</th>\n",
       "      <th>总出险数量</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>-2.5601</td>\n",
       "      <td>680.1081</td>\n",
       "      <td>2897297.4500</td>\n",
       "      <td>1236.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>-2.6360</td>\n",
       "      <td>704.3563</td>\n",
       "      <td>4185235.5200</td>\n",
       "      <td>1282.0000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        V1       V2        总理赔金额     总出险数量\n",
       "36 -2.5601 680.1081 2897297.4500 1236.0000\n",
       "37 -2.6360 704.3563 4185235.5200 1282.0000"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 真实值如下\n",
    "df_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a089785-db0c-4a0e-beef-f9e2d2c50a43",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 最后预测\n",
    "# 选择lag=4拟合模型\n",
    "model_fitted = model.fit(4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
